A numerical framework for linear stability analysis of two-phase stratified pipe flows

TL;DR

This paper introduces a novel numerical framework for linear stability analysis of two-phase stratified pipe flows, considering three-dimensional disturbances including interface effects, verified through Cartesian coordinate solutions, applicable to gas-liquid and liquid-liquid flows.

Contribution

It presents the first comprehensive 3D stability analysis method for two-phase pipe flows using bipolar coordinates, enabling detailed interface and shear stress investigations.

Findings

Validated bipolar coordinate method against Cartesian solutions.

Analyzed stability of gas-liquid and liquid-liquid flows.

Demonstrated detailed flow behavior near the interface and wall.

Abstract

A numerical framework for rigorous linear stability analysis of two-phase stratified flows of two immiscible fluids in horizontal circular pipes is presented. For the first time, three-dimensional disturbances, including those at the interface between two fluids, are considered. The proposed numerical framework is based on a finite-volume method and allows solving the problem numerically in bipolar cylindrical coordinates. In these coordinates, both the pipe wall and the unperturbed interface (of a constant curvature, e.g., plane interface, as considered in this work) coincide with the coordinate surfaces. Thereby, the no-slip as well as the interfacial boundary conditions can be imposed easily. It also enables investigation of the local behavior of the flow field and shear stresses in the vicinity of the triple points, where the interface contacts the pipe wall. The results obtained in the bipolar coordinates are verified by an independent numerical solution based on the problem formulation in Cartesian coordinates, where the pipe wall is treated by the immersed boundary method. Two representative examples of gas-liquid and liquid-liquid flows are included to demonstrate the applicability of the proposed numerical technique for analyzing the flow stability.

Tables

Table 1. Table 1: Physical properties of the two-phase pipe flows considered.

	${\mu }_1$	${\mu }_2$	${\mu }_{12}$	${\rho }_1$	${\rho }_2$	${\rho }_{12}$	$\sigma$	$D$
	$\left[\frac{kg}{m\cdot s}\right]$	$\left[\frac{kg}{m\cdot s}\right]$		$\left[\frac{kg}{m^3}\right]$	$\left[\frac{kg}{m^3}\right]$		$\left[\frac{N}{m}\right]$	$\left[m\right]$
oil-water	0.001	0.0097	0.103	1000	835	1.197	0.03	0.02
air-water	0.001	1.818$\times {10}^{-5}$	55.556	1000	1	1000	0.072	0.014

Table 2. Table 2: Base flow convergence. Liquid-liquid ( m = 0.103 𝑚 0.103 m=0.103 ) flow in a circular pipe of diameter D = 0.02 𝐷 0.02 D=0.02 m.

		$U_{1S}/U_{2S}$
$N_{\xi }$	$N_{\varphi }$	$h=0.202$		$h=0.601$
$\left(N_x\right)$	$\left(N_y\right)$	Bipolar FVM	IBM	Bipolar FVM	IBM
50	50	0.3813	0.3828	6.2936	6.3054
100	100	0.3801	0.3808	6.2821	6.2815
200	200	0.3797	0.3800	6.2788	6.2787
300	300	0.3796	0.3798	6.2782	6.2782
400	400	0.3796	0.3797	6.2780	6.2780
Analytical		0.3801		6.2864

Table 3. Table 3: Critical flow rates and the critical eigenvalue ( α c ​ r = 5.6 subscript 𝛼 𝑐 𝑟 5.6 \alpha_{cr}=5.6 ) for h = 0.202 ℎ 0.202 h=0.202 for liquid-liquid ( m = 0.103 𝑚 0.103 m=0.103 , r = 1.198 𝑟 1.198 r=1.198 ) flow in a circular pipe of diameter D = 0.02 𝐷 0.02 D=0.02 m.

$N_{\xi }$	$N_{\varphi }$	Bipolar FVM			IBM
$\left(N_x\right)$	$\left(N_y\right)$	$U_{1S}$	$U_{2S}$	${\lambda }_I$	$U_{1S}$	$U_{2S}$	${\lambda }_I$
50	50	0.0653	0.1720	-10.9131	0.0759	0.2000	-10.5352
100	100	0.0664	0.1751	-10.8338	0.0699	0.1841	-10.7157
200	200	0.0666	0.1755	-10.8297	0.0678	0.1786	-10.7841
300	300	0.0667	0.1756	-10.8293	0.0673	0.1773	-10.8021
400	400	0.0667	0.1756	-10.8293	0.0672	0.1769	-10.8085

Table 4. Table 4: Holdup and critical eigenvalue for different ξ max subscript 𝜉 max \xi_{\text{max}} . Liquid-liquid ( m = 0.103 𝑚 0.103 m=0.103 , r = 1.198 𝑟 1.198 r=1.198 ) flow in a circular pipe of diameter D = 0.02 𝐷 0.02 D=0.02 m, U 1 ​ S = 0.0666 subscript 𝑈 1 𝑆 0.0666 U_{1S}=0.0666 , U 2 ​ S = 0.1755 subscript 𝑈 2 𝑆 0.1755 U_{2S}=0.1755 , α = 5.6 𝛼 5.6 \alpha=5.6 , N ξ = 200 subscript 𝑁 𝜉 200 N_{\xi}=200 , N ϕ = 200 subscript 𝑁 italic-ϕ 200 N_{\phi}=200 .

${\xi }_{\text{max}}$	$h$	${\lambda }_R$	${\lambda }_I$
$\pi$	0.2017	0.0201	-10.7683
$2\pi$	0.2021	-0.7070 $\times {10}^{-4}$	-10.8294
$3\pi$	0.2021	0.5632 $\times {10}^{-10}$	-10.8297
$4\pi$	0.2021	0.8255 $\times {10}^{-4}$	-10.8295
$5\pi$	0.2021	0.000144	-10.8293

Table 5. Table 5: Long-wave critical superficial velocities and critical eigenvalue for h = 0.6 ℎ 0.6 h=0.6 . Air-water ( m = 55 𝑚 55 m=55 , r = 1000 𝑟 1000 r=1000 ) flow in a circular pipe of diameter D = 0.014 𝐷 0.014 D=0.014 m ( N ξ = N ϕ = 200 subscript 𝑁 𝜉 subscript 𝑁 italic-ϕ 200 N_{\xi}=N_{\phi}=200 ).

$\alpha$	$U_{1S}$	$U_{2S}$	${\lambda }_I$	$c$
0.01	0.0488	0.3247	-0.010863	1.0862
0.001	0.0481	0.3197	-0.001098	1.0975
0.0001	0.0481	0.3197	-0.000110	1.0976
0.00001	0.0481	0.3197	-0.1097$\times {10}^{-4}$	1.0976

Table 6. Table 6: Critical superficial velocities and the critical eigenvalue ( α c ​ r = 0.001 subscript 𝛼 𝑐 𝑟 0.001 \alpha_{cr}=0.001 ) for h = 0.6 ℎ 0.6 h=0.6 . Air-water ( m = 55 𝑚 55 m=55 , r = 1000 𝑟 1000 r=1000 ) flow in a circular pipe of diameter D = 0.014 𝐷 0.014 D=0.014 m.

$N_{\xi }$	$N_{\varphi }$	Bipolar FVM			IBM
$\left(N_x\right)$	$\left(N_y\right)$	$U_{1S}$	$U_{2S}$	${\lambda }_I$	$U_{1S}$	$U_{2S}$	${\lambda }_I$
50	50	0.0476	0.3168	-0.001085	0.0504	0.3351	-0.001062
100	100	0.0480	0.3191	-0.001098	0.0503	0.3348	-0.001064
200	200	0.0481	0.3197	-0.001098	0.0504	0.3352	-0.001063
300	300	0.0481	0.3198	-0.001097	0.0504	0.3352	-0.001062
400	400	0.0481	0.3199	-0.001098	0.0504	0.3353	-0.001062

Equations

$$\nabla\cdot{\mbox{\boldmath$u$}}^{(k)}=0,$$

$$\frac{\partial{\mbox{\boldmath$u$}}^{(k)}}{\partial t}+\left({\mbox{\boldmath$u$}}^{(k)}\cdot\nabla\right){\mbox{\boldmath$u$}}^{(k)}=-\frac{\rho_{1}}{\rho_{12}\rho_{k}}\nabla p^{(k)}+\frac{1}{\mbox{Re}}\frac{\rho_{1}}{\rho_{12}\rho_{k}}\frac{\mu_{12}\mu_{k}}{\mu_{1}}\Delta{\mbox{\boldmath$u$}}^{(k)}+\frac{1}{\mbox{Fr}}{\mbox{\boldmath$e_{g}$}},$$

$$x$$

$$y$$

$$\xi$$

$$\varphi$$

$$H_{\xi }=H_{\varphi }=\frac{1}{2}\frac{\sin {\varphi }_0}{\cosh \xi -\cos \varphi }.$$

$$\frac{1}{H_{\xi }}\frac{\partial u_{\xi }^{(k)}}{\partial \xi }+\frac{1}{H_{\varphi }}\frac{\partial u_{\varphi }^{(k)}}{\partial \varphi }-2\frac{\sinh \xi }{\sin {\varphi }_0}{u}_{\xi }^{(k)}-2\frac{\sin \varphi }{\sin {\varphi }_0}{u}_{\varphi }^{(k)}+\frac{\partial u_z}{\partial z}=0,$$

$$\displaystyle\begin{aligned} &\frac{\partial u_{\xi}^{(k)}}{\partial t}+\frac{u_{\xi}^{(k)}}{H_{\xi}}\frac{\partial u_{\xi}^{(k)}}{\partial\xi}+\frac{u_{\phi}^{(k)}}{H_{\phi}}\frac{\partial u_{\xi}^{(k)}}{\partial\phi}+u_{z}^{(k)}\frac{\partial u_{\xi}^{(k)}}{\partial z}-2\frac{\sin\phi}{\sin\phi_{0}}u_{\xi}^{(k)}u_{\phi}^{(k)}+2\frac{\sinh\xi}{\sin\phi_{0}}\left(u_{\phi}^{(k)}\right)^{2}\\ &=-\frac{1}{\rho_{k}}\frac{1}{H_{\xi}}\frac{\partial p^{(k)}}{\partial\xi}+\frac{1}{\mbox{Re}}\frac{\rho_{1}}{\rho_{12}\rho_{k}}\frac{\mu_{12}\mu_{k}}{\mu_{1}}\Biggl{[}\frac{1}{H_{\xi}H_{\phi}}\frac{\partial}{\partial\xi}\biggl{(}\frac{H_{\phi}}{H_{\xi}}\frac{\partial u_{\xi}^{(k)}}{\partial\xi}\biggr{)}+\frac{1}{H_{\xi}H_{\phi}}\frac{\partial}{\partial\phi}\biggl{(}\frac{H_{\xi}}{H_{\phi}}\frac{\partial u_{\xi}^{(k)}}{\partial\phi}\biggr{)}\\ &+\frac{\partial^{2}u_{\xi}^{(k)}}{\partial z^{2}}-4\frac{\sin\phi}{\sin\phi_{0}}\frac{1}{H_{\xi}}\frac{\partial u_{\phi}^{(k)}}{\partial\xi}+4\frac{\sinh\xi}{\sin\phi_{0}}\frac{1}{H_{\phi}}\frac{\partial u_{\phi}^{(k)}}{\partial\phi}+4\frac{\cos^{2}\phi-\cosh^{2}\xi}{\sin^{2}\phi_{0}}u_{\xi}^{(k)}\Biggr{]}\\ &+\frac{1}{\mbox{Fr}}\frac{\sinh\xi\sin\phi}{\bigl{(}\cosh\xi-\cos\phi\bigr{)}},\end{aligned}$$

$$\displaystyle\begin{aligned} \frac{\partial u_{\phi}^{(k)}}{\partial t}&+\frac{u_{\xi}^{(k)}}{H_{\xi}}\frac{\partial u_{\phi}^{(k)}}{\partial\xi}+\frac{u_{\phi}^{(k)}}{H_{\phi}}\frac{\partial u_{\phi}^{(k)}}{\partial\phi}+u_{z}^{(k)}\frac{\partial u_{\phi}^{(k)}}{\partial z}-2\frac{\sinh\xi}{\sin\phi_{0}}u_{\xi}^{(k)}u_{\phi}^{(k)}+2\frac{\sin\phi}{\sin\phi_{0}}\left(u_{\xi}^{(k)}\right)^{2}\\ =&-\frac{1}{\rho_{k}}\frac{1}{H_{\phi}}\frac{\partial p^{(k)}}{\partial\phi}+\frac{1}{\mbox{Re}}\frac{\rho_{1}}{\rho_{12}\rho_{k}}\frac{\mu_{12}\mu_{k}}{\mu_{1}}\Biggl{[}\frac{1}{H_{\xi}H_{\phi}}\frac{\partial}{\partial\xi}\biggl{(}\frac{H_{\phi}}{H_{\xi}}\frac{\partial u_{\phi}^{(k)}}{\partial\xi}\biggr{)}\\ &+\frac{1}{H_{\xi}H_{\phi}}\frac{\partial}{\partial\phi}\biggl{(}\frac{H_{\xi}}{H_{\phi}}\frac{\partial u_{\phi}^{(k)}}{\partial\phi}\biggr{)}+\frac{\partial^{2}u_{\phi}^{(k)}}{\partial z^{2}}+4\frac{\sin\phi}{\sin\phi_{0}}\frac{1}{H_{\xi}}\frac{\partial u_{\xi}^{(k)}}{\partial\xi}-4\frac{\sinh\xi}{\sin\phi_{0}}\frac{1}{H_{\phi}}\frac{\partial u_{\xi}^{(k)}}{\partial\phi}\\ &+4\frac{\cos^{2}\phi-\cosh^{2}\xi}{\sin^{2}\phi_{0}}u_{\phi}^{(k)}\Biggr{]}+\frac{1}{\mbox{Fr}}\frac{\bigl{(}1-\cosh\xi\cos\phi\bigr{)}}{\bigl{(}\cosh\xi-\cos\phi\bigr{)}},\end{aligned}$$

$$\displaystyle\begin{aligned} \frac{\partial u_{z}^{(k)}}{\partial t}&+\frac{u_{\xi}^{(k)}}{H_{\xi}}\frac{\partial u_{z}^{(k)}}{\partial\xi}+\frac{u_{\phi}^{(k)}}{H_{\phi}}\frac{\partial u_{z}^{(k)}}{\partial\phi}+u_{z}^{(k)}\frac{\partial u_{z}^{(k)}}{\partial z}=-\frac{1}{\rho_{k}}\frac{\partial p^{(k)}}{\partial z}\\ &+\frac{1}{\mbox{Re}}\frac{\rho_{1}}{\rho_{12}\rho_{k}}\frac{\mu_{12}\mu_{k}}{\mu_{1}}\Biggl{[}\frac{1}{H_{\xi}H_{\phi}}\frac{\partial}{\partial\xi}\left(\frac{H_{\phi}}{H_{\xi}}\frac{\partial u_{z}^{(k)}}{\partial\xi}\right)+\frac{1}{H_{\xi}H_{\phi}}\frac{\partial}{\partial\phi}\left(\frac{H_{\xi}}{H_{\phi}}\frac{\partial u_{z}^{(k)}}{\partial\phi}\right)\\ &+\frac{\partial^{2}u_{z}^{(k)}}{\partial z^{2}}\Biggr{]}.\end{aligned}$$

$${\mbox{\boldmath$u$}}^{(1)}(\phi=\phi_{0}+\pi)={\mbox{\boldmath$u$}}^{(2)}(\phi=\phi_{0})=0$$

$${\varphi }_{\eta }=\pi +\eta \left(\xi ,z\right).$$

$$u_{\varphi }^{(k)}(\varphi ={\varphi }_{\eta })=H_{\varphi }\frac{\partial \eta }{\partial t}+u_{\xi }^{(k)}(\varphi ={\varphi }_{\eta })\frac{H_{\varphi }}{H_{\xi }}\frac{\partial \eta }{\partial \xi }+u_z^{(k)}(\varphi ={\varphi }_{\eta })H_{\varphi }\frac{\partial \eta }{\partial z}.$$

$${\mbox{\boldmath$u$}}^{(1)}(\phi=\phi_{\eta})={\mbox{\boldmath$u$}}^{(2)}(\phi=\phi_{\eta}).$$

$${\tau }_{\xi \varphi }^{(1)}(\varphi ={\varphi }_{\eta })$$

$${\tau }_{\varphi z}^{(1)}(\varphi ={\varphi }_{\eta })$$

$$\displaystyle\biggl{\llbracket}-p+\sigma_{n}\biggr{\rrbracket}_{\phi=\phi_{\eta}}$$

$$\frac{1}{H_{\varphi }{H}_{\xi }}\left[\frac{\partial }{\partial \xi }\left(\frac{H_{\varphi }}{H_{\xi }}\frac{\partial U_z^{(k)}}{\partial \xi }\right)+\frac{\partial }{\partial \varphi }\left(\frac{H_{\xi }}{H_{\varphi }}\frac{\partial U_z^{(k)}}{\partial \varphi }\right)\right]$$

$$-\frac{\partial P^{(k)}}{\partial \xi }+\frac{{\rho }_k}{\text{mbox}Fr}{H}_{\xi }\frac{\sinh \xi \sin \varphi }{\cosh \xi -\cos \varphi }$$

$$\displaystyle-\frac{\partial P^{(k)}}{\partial\phi}+\frac{\rho_{k}}{\mbox{Fr}}H_{\phi}\frac{\bigl{(}1-\cosh\xi\cos\phi\bigr{)}}{\cosh\xi-\cos\phi}$$

$$U_z^{(1)}\left(\varphi ={\varphi }_0+\pi \right)=0,U_z^{(2)}\left(\varphi ={\varphi }_0\right)=0,$$

$$U_z^{(1)}(\varphi =\pi )=U_z^{(2)}(\varphi =\pi ),$$

$$\mu_{12}\bigg{(}\frac{\partial U_{z}^{(1)}}{\partial\phi}\bigg{)}_{\phi=\pi}=\bigg{(}\frac{\partial U_{z}^{(2)}}{\partial\phi}\bigg{)}_{\phi=\pi}.$$

$$P^{(1)}(\varphi =\pi )=P^{(2)}(\varphi =\pi ).$$

$${\int }_{\pi }^{{\varphi }_0+\pi }{\int }_{-\infty }^{+\infty }{U}_z^{(1)}{H}_{\xi }{H}_{\varphi }d\xi d\varphi$$

$${\int }_{{\varphi }_0}^{\pi }{\int }_{-\infty }^{+\infty }{U}_z^{(2)}{H}_{\xi }{H}_{\varphi }d\xi d\varphi$$

$$u^{(k)}$$

$$p^{(k)}$$

$$\begin{pmatrix}{\mbox{\boldmath$u^{\prime(k)}$}}\\ p^{\prime(k)}\\ \eta\end{pmatrix}=\begin{pmatrix}{\mbox{\boldmath$\tilde{u}^{(k)}$}}(\xi,\phi)=\left[\tilde{u}_{\xi}^{(k)}(\xi,\phi),\tilde{u}_{\phi}^{(k)}(\xi,\phi),\tilde{u}_{z}^{(k)}(\xi,\phi)\right]\\ \tilde{p}^{(k)}(\xi,\phi)\\ \tilde{\eta}(\xi)\end{pmatrix}e^{\left(i\alpha z+\lambda t\right)},$$

$$\displaystyle\begin{aligned} f({\mbox{\boldmath$r$}})&=f({\mbox{\boldmath$r_{0}$}})+\nabla f\biggr{|}_{{\mbox{\boldmath$r_{0}$}}}\cdot\left({\mbox{\boldmath$r$}}-{\mbox{\boldmath$r_{0}$}}\right)+\mathcal{O}(\left({\mbox{\boldmath$r$}}-{\mbox{\boldmath$r_{0}$}}\right)^{2})=f({\mbox{\boldmath$r_{0}$}})+\biggl{(}\frac{1}{H_{\xi}}\frac{\partial f}{\partial\xi}\biggr{|}_{{\mbox{\boldmath$r_{0}$}}}{\mbox{\boldmath$e_{\xi}$}}\\ &+\frac{1}{H_{\phi}}\frac{\partial f}{\partial\phi}\biggr{|}_{{\mbox{\boldmath$r_{0}$}}}{\mbox{\boldmath$e_{\phi}$}}+\frac{\partial f}{\partial z}\biggr{|}_{{\mbox{\boldmath$r_{0}$}}}{\mbox{\boldmath$e_{z}$}}\biggr{)}\cdot\bigg{(}0{\mbox{\boldmath$e_{\xi}$}}+\eta\left(\xi,z,t\right)H_{\phi}{\mbox{\boldmath$e_{\phi}$}}+0{\mbox{\boldmath$e_{z}$}}\bigg{)}+\mathcal{O}(\eta^{2})\\ &=f\bigl{(}\xi,\pi,z\bigr{)}+\frac{1}{H_{\phi}}\frac{\partial f}{\partial\phi}\biggr{|}_{\phi=\pi}H_{\phi}\eta\left(\xi,z,t\right)+\mathcal{O}(\eta^{2}),\end{aligned}$$

$$\frac{1}{H_{\xi }}\frac{\partial {\tilde{u}}_{\xi }^{(k)}}{\partial \xi }+\frac{1}{H_{\varphi }}\frac{\partial {\tilde{u}}_{\varphi }^{(k)}}{\partial \varphi }-2\frac{\sinh \xi }{\sin {\varphi }_0}{\tilde{u}}_{\xi }^{(k)}-2\frac{\sin \varphi }{\sin {\varphi }_0}{\tilde{u}}_{\varphi }^{(k)}+i\alpha {\tilde{u}}_z=0,$$

$$\displaystyle\begin{aligned} \lambda\tilde{u}_{\xi}^{(k)}&+i\alpha U_{z}^{(k)}\tilde{u}_{\xi}^{(k)}=-\frac{\rho_{1}}{\rho_{12}\rho_{k}}\frac{1}{H_{\xi}}\frac{\partial\tilde{p}^{(k)}}{\partial\xi}\\ &+\frac{1}{\mbox{Re}}\frac{\rho_{1}}{\rho_{12}\rho_{k}}\frac{\mu_{12}\mu_{k}}{\mu_{1}}\Biggl{[}\frac{1}{H_{\xi}H_{\phi}}\frac{\partial}{\partial\xi}\biggl{(}\frac{H_{\phi}}{H_{\xi}}\frac{\partial\tilde{u}_{\xi}^{(k)}}{\partial\xi}\biggr{)}+\frac{1}{H_{\xi}H_{\phi}}\frac{\partial}{\partial\phi}\biggl{(}\frac{H_{\xi}}{H_{\phi}}\frac{\partial\tilde{u}_{\xi}^{(k)}}{\partial\phi}\biggr{)}-\alpha^{2}\tilde{u}_{\xi}^{(k)}\\ &-4\frac{\sin\phi}{\sin\phi_{0}}\frac{1}{H_{\xi}}\frac{\partial\tilde{u}_{\phi}^{(k)}}{\partial\xi}+4\frac{\sinh\xi}{\sin\phi_{0}}\frac{1}{H_{\phi}}\frac{\partial\tilde{u}_{\phi}^{(k)}}{\partial\phi}+4\frac{\cos^{2}\phi-\cosh^{2}\xi}{\sin^{2}\phi_{0}}\tilde{u}_{\xi}^{(k)}\Biggr{]},\end{aligned}$$

$$\displaystyle\begin{aligned} \lambda\tilde{u}_{\phi}^{(k)}&+i\alpha U_{z}^{(k)}\tilde{u}_{\phi}^{(k)}=-\frac{\rho_{1}}{\rho_{12}\rho_{k}}\frac{1}{H_{\xi}}\frac{\partial\tilde{p}^{(k)}}{\partial\phi}\\ &+\frac{1}{\mbox{Re}}\frac{\rho_{1}}{\rho_{12}\rho_{k}}\frac{\mu_{12}\mu_{k}}{\mu_{1}}\Biggl{[}\frac{1}{H_{\xi}H_{\phi}}\frac{\partial}{\partial\xi}\biggl{(}\frac{H_{\phi}}{H_{\xi}}\frac{\partial\tilde{u}_{\phi}^{(k)}}{\partial\xi}\biggr{)}+\frac{1}{H_{\xi}H_{\phi}}\frac{\partial}{\partial\phi}\biggl{(}\frac{H_{\xi}}{H_{\phi}}\frac{\partial\tilde{u}_{\phi}^{(k)}}{\partial\phi}\biggr{)}-\alpha^{2}\tilde{u}_{\phi}^{(k)}\\ &+4\frac{\sin\phi}{\sin\phi_{0}}\frac{1}{H_{\xi}}\frac{\partial\tilde{u}_{\xi}^{(k)}}{\partial\xi}-4\frac{\sinh\xi}{\sin\phi_{0}}\frac{1}{H_{\phi}}\frac{\partial\tilde{u}_{\xi}^{(k)}}{\partial\phi}+4\frac{\cos^{2}\phi-\cosh^{2}\xi}{\sin^{2}\phi_{0}}\tilde{u}_{\phi}^{(k)}\Biggr{]},\end{aligned}$$

$$\displaystyle\begin{aligned} \lambda\tilde{u}_{z}^{(k)}&+\frac{\tilde{u}_{\xi}^{(k)}}{H_{\xi}}\frac{\partial U_{z}^{(k)}}{\partial\xi}+\frac{\tilde{u}_{\phi}^{(k)}}{H_{\phi}}\frac{\partial U_{z}^{(k)}}{\partial\phi}+i\alpha U_{z}^{(k)}\tilde{u}_{z}^{(k)}=-i\alpha\frac{\rho_{1}}{\rho_{12}\rho_{k}}\tilde{p}^{(k)}\\ &+\frac{1}{\mbox{Re}}\frac{\rho_{1}}{\rho_{12}\rho_{k}}\frac{\mu_{12}\mu_{k}}{\mu_{1}}\Biggl{[}\frac{1}{H_{\xi}H_{\phi}}\frac{\partial}{\partial\xi}\left(\frac{H_{\phi}}{H_{\xi}}\frac{\partial\tilde{u}_{z}^{(k)}}{\partial\xi}\right)+\frac{1}{H_{\xi}H_{\phi}}\frac{\partial}{\partial\phi}\left(\frac{H_{\xi}}{H_{\phi}}\frac{\partial\tilde{u}_{z}^{(k)}}{\partial\phi}\right)\\ &-\alpha^{2}\tilde{u}_{z}^{(k)}\Biggr{]}.\end{aligned}$$

$${\mbox{\boldmath$\tilde{u}$}}^{(1)}(\phi=\phi_{0}+\pi)={\mbox{\boldmath$\tilde{u}$}}^{(2)}(\phi=\phi_{0})=0.$$

Full text

A numerical framework for linear stability analysis of two-phase stratified pipe flows

Ilya Barmak111Corresponding author. E-mail address:

ilyab@tauex.tau.ac.il

Alexander Gelfgat

Neima Brauner

Abstract

A numerical framework for rigorous linear stability analysis of two-phase stratified flows of two immiscible fluids in horizontal circular pipes is presented. For the first time, three-dimensional disturbances, including those at the interface between two fluids, are considered. The proposed numerical framework is based on a finite–volume method and allows solving the problem numerically in bipolar cylindrical coordinates. In these coordinates, both the pipe wall and the unperturbed interface (of a constant curvature, e.g., plane interface, as considered in this work) coincide with the coordinate surfaces. Thereby, the no-slip as well as the interfacial boundary conditions can be imposed easily. It also enables investigation of the local behavior of the flow field and shear stresses in the vicinity of the triple points, where the interface contacts the pipe wall. The results obtained in the bipolar coordinates are verified by an independent numerical solution based on the problem formulation in Cartesian coordinates, where the pipe wall is treated by the immersed boundary method. Two representative examples of gas-liquid and liquid-liquid flows are included to demonstrate the applicability of the proposed numerical technique for analyzing the flow stability.

1 Introduction

In many industrial processes involving two-phase liquid-liquid and gas-liquid pressure-driven flows (e.g., gas-condensate pipeline transport), the desirable flow configuration is the stratified smooth flow pattern, in which heavy and light fluids form two continuous layers, lower and upper, respectively. The stratified smooth flow pattern, however, is feasible only for certain operational conditions for which the flow is stable and the fluid-fluid interface remains smooth. The critical conditions corresponding to the onset of interfacial perturbations have been explored experimentally (e.g., Charles and Lilleleht,, 1965; Yu and Sparrow,, 1969; Kao and Park,, 1972; Barnea et al.,, 1980; Kokal and Stanislav,, 1989) and theoretically. To avoid complexity of rigorous stability analysis in circular pipe geometry, many theoretical studies were based on simplified mechanistic models, mainly in the framework of the Two-Fluid (TF) model (e.g., Andritsos et al.,, 1989; Brauner and Maron,, 1991; Barnea and Taitel,, 1993; Ullmann and Brauner,, 2006; Kushnir et al.,, 2017).

In other studies, the comprehensive linear stability analysis was applied to two-phase flow between two infinite plates. This simplified configuration allows one to account for all possible infinitesimal perturbations (e.g., Yih,, 1967; Charru and Fabre,, 1994; Ó Náraigh et al.,, 2014; Kaffel and Riaz,, 2015; Barmak et al.,, 2016, 2019). This approach yielded better understanding of mechanisms involved in the onset of instability of stratified flows. However, the simplification of geometry allows for only qualitative comparison of the stability characteristics with experimental data obtained in pipe flow.

The realistic cross-sectional geometry of a rectangular duct was used recently in stability studies of horizontal and inclined two-phase stratified flows (Gelfgat and Brauner,, 2020; Gelfgat et al.,, 2021). A qualitative agreement of stability characteristics with the flows in the two-plate geometry (TP) was established. At the same time, the quantitative differences were found to be substantial. The effect of circular pipe walls on the onset of instability in two-phase stratified flows has never been studied, and it is still unknown whether the effect is at least qualitatively similar to that for either the TP or square duct geometry. An additional complication found in two-phase pipe flows is that in the vicinity of triple points (Goldstein and Eyal,, 2021), i.e., the contact points of the fluid–fluid interface with the pipe wall, the shear stresses may differ strongly from those in the flow bulk (e.g., Brauner et al.,, 1996; Goldstein et al.,, 2021).

It is well-known that the single-phase Poiseuille flow in a circular pipe is linearly stable for all Reynolds numbers (e.g. Davey and Drazin,, 1969), while Poiseuille flow between two plates becomes unstable at Re=5772 (Orszag,, 1971). Since the transition to turbulence takes place at the Reynolds numbers slightly above 2000, the linear stability results are irrelevant for predicting transitions in these flows. Contrarily, the linear stability analysis describes the primary instability in the two-phase flows correctly, as was recently shown by comparison of experimental and numerical results for square ducts in (Nezihovski et al.,, 2022). Therefore, one can expect that linear stability analysis conducted for a two-phase flow in a circular pipe will also describe the primary transition correctly. To the best of our knowledge, no attempt has been made to carry out a rigorous stability analysis of stratified two-phase flows in circular pipes, which is the main objective of this study.

In the present study, we formulate and solve the linear stability problem in bipolar cylindrical coordinates (Korn and Korn,, 2000), where both a constant curvature interface and the pipe wall coincide with the coordinate surfaces. In these coordinates, the boundary conditions can be imposed relatively easily and there is an analytical solution available for the steady fully developed flow (e.g., Goldstein et al.,, 2015). We also apply an alternative approach, treating the pipe wall by the immersed boundary method (IBM) in Cartesian coordinates, in which the plane interface coincides with the plane of constant vertical coordinate. Since there are no independent results on stability of two-phase stratified pipe flow to compare with, the results obtained by two approaches cross verify each other.

In the following, we formulate the problem in the vector form and in the bipolar coordinates. We present in detail the finite-volume method in bipolar coordinates, which appears to be non-trivial. For two representative examples of gas-liquid and liquid-liquid flows, we perform a comprehensive stability analysis and discuss the stability characteristics.

2 Problem formulation

Two-phase stratified flow in a horizontal circular pipe, which is driven by pressure gradient imposed in the direction of the pipe axis, is a particular case of a general two-layer flow of two immiscible fluids in a channel of constant cross section. For such a flow, we formulate governing equations both in a vector form and in specially tailored curvilinear orthogonal coordinates that fit the boundaries of the problem, i.e., both the channel walls and the unperturbed interface coincide with coordinate lines (Appendix A).

2.1 Governing equations in vector form

The two-phase flow is described by velocity {\mbox{\boldmathu}}^{(k)} and pressure $p^{(k)}$ fields that satisfy the dimensionless continuity and momentum equations defined in each layer $k=1,2$ (1 - heavy fluid, 2 - light fluid):

[TABLE]

[TABLE]

where the velocity scale is the mixture velocity $\displaystyle U_{m}=U_{1S}+U_{2S}$. $U_{1S}$ and $U_{2S}$ represent the superficial velocities of each phase, i.e., the flow rate of each phase divided by the cross-sectional area of the channel. The time and pressure are scaled by $\displaystyle D/U_{m}$ ($D$ - characteristic channel size, e.g., pipe diameter) and $\rho_{2}U_{m}^{2}$, respectively. The dimensionless governing parameters are the density and viscosity ratios, $\displaystyle\rho_{12}=\rho_{1}/\rho_{2}$ and $\displaystyle\mu_{12}=\mu_{1}/\mu_{2}$, respectively, the Reynolds number $\displaystyle\mbox{Re}=\rho_{2}DU_{m}/\mu_{2}$ and the Froude number $\displaystyle\mbox{Fr}=U_{m}^{2}/(gD)$, where $g$ is the gravitational acceleration. $e_{g}$ is a unit vector in the direction of gravity.

2.2 Flow geometry

For the considered two-phase stratified flow in a circular pipe of diameter $D$, the convenient coordinates are the bipolar cylindrical coordinates $(\xi,\phi,z)$ (e.g., Korn and Korn,, 2000; Goldstein et al.,, 2015), in which two lines of foci align with the contact lines of two phases and pipe wall, while the z-axis coincides with the pipe axis. In the pipe cross section (Fig. 1), $\xi$ is defined as the natural logarithm of the ratio of distances from a point to the foci F1 and F2 and $\phi$ as the angle formed by the two foci and the point. For example, for a point P on the pipe wall, $\displaystyle\xi_{P}=\ln(|\text{PF}_{2}|/|\text{PF}_{1}|)$ and $\displaystyle\phi_{P}=\phi_{0}=\angle\text{F}_{1}\text{PF}_{2}$. Thus, in the bipolar coordinates, the pipe wall corresponds to isolines of coordinate $\phi$. The upper section of the pipe wall bounds the light phase (shown in red in Fig. 1) and is represented by $\phi_{0}$, while the bottom of the pipe, bounding the heavy phase (shown in blue in Fig. 1), is represented by $\phi_{0}+\pi$. The interface that separates two fluids is a surface defined as $\phi_{\eta}=\phi^{*}+\eta\left(\xi,z\right)$.

In this work, we consider gravity-dominated systems with Eötvös number, $\displaystyle\mbox{Eo}=(\rho_{1}-\rho_{2})gD^{2}/\sigma\gg 1$, where $\sigma$ is the surface tension coefficient. The smooth stratified flow in such systems is characterized by a plane interface (Gorelik and Brauner,, 1999), so that $\phi^{*}=\pi$ (as shown in Fig. 1). The heavy fluid occupies an infinite strip of \displaystyle\bigl{(}-\infty<\xi<+\infty,\pi<\phi<\phi_{0}+\pi\bigr{)} with the cross-sectional area of \displaystyle A_{1}=0.25D^{2}\big{(}\phi_{0}-0.5\sin 2\phi_{0}\big{)}, while the light fluid – an infinite strip of \displaystyle\bigl{(}-\infty<\xi<+\infty,\phi_{0}<\phi<\pi\bigr{)} with the area of $\displaystyle A_{2}=A-A_{1}$, where $\displaystyle A=\pi D^{2}/4$ is the cross-sectional area of the pipe. Then, one of the dimensionless parameters defining two-phase stratified flow is the heavy phase holdup, defined as $\displaystyle h=A_{1}/A$. In the following, we normalize all the length scales by the pipe diameter, $D$.

In Cartesian coordinates, the pipe center is located at $x,y=(0,0)$, and the foci of the bipolar coordinate system are located in the triple points, $\displaystyle\text{F}_{1}(-0.5\sin\phi_{0},-0.5\cos\phi_{0})$ and $\displaystyle\text{F}_{2}(0.5\sin\phi_{0},-0.5\cos\phi_{0})$. The relations of the Cartesian coordinates to the bipolar coordinates are

[TABLE]

The reciprocal relations from the Cartesian to the bipolar coordinates are

[TABLE]

The Lamé coefficients (scale factors) for the bipolar coordinates are

[TABLE]

It is worth noting that the cylindrical coordinate system, being a natural choice for computation of a single-phase flow in a circular pipe, is not convenient for the modeling of a two-phase stratified pipe flow. In addition to possible cross flows across the cylindrical axis, in two-phase flows, the fluid—fluid interface can intersect the axis, which is problematic due to its singularity in numerical formulation. The second problem in using cylindrical coordinates is the treatment of the interface. In case it is assumed to be plane in the unperturbed state (valid for gravity-dominated systems, e.g., Gorelik and Brauner,, 1999), the interface depends on two cross-sectional coordinates and the formulation and implementation of the non-trivial interfacial boundary conditions becomes intricate. Moreover, as we will show below, the interface can be located at any height, depending on the flow rates of two fluids. In the following, we provide a formulation in bipolar cylindrical coordinates that appear to be the most suitable.

2.3 Governing equations in bipolar coordinates

The governing equations for three-dimensional time-dependent flow in the bipolar cylindrical coordinates can be obtained from the corresponding equations in general orthogonal coordinates, Eqs. A.2 and A.3 (see Appendix A), by introducing the Lamé coefficients (Eq. 5):

[TABLE]

[TABLE]

Note that $\displaystyle H_{\xi}d\xi$ and $\displaystyle H_{\phi}d\phi$ have the meaning of physical lengths so that they are important for the purpose of the finite–volume integration on the computational grid and should not be canceled out (see sec. 5).

2.4 Boundary conditions

The governing equations for the velocity field are coupled with the no-slip conditions at the pipe wall, including at the triple points ($\xi\to\pm\infty$)

[TABLE]

and the boundary conditions at the interface between the two fluids. Here we present only the final expressions for the interfacial boundary conditions, while the detailed description can be found in Appendix B. The interface, which is considered plane in an unperturbed state, is described as

[TABLE]

The absence of mass transfer across the interface is accounted by the kinematic boundary condition

[TABLE]

The velocity field is continuous across the interface, therefore

[TABLE]

The interfacial boundary conditions require continuity of the tangential components of the viscous stress tensor, i.e., $\displaystyle\tau_{\xi\phi}$ in the plane $\displaystyle(\xi,\phi)$ and $\displaystyle\tau_{\phi z}$ in the plane $\displaystyle(\phi,z)$, while the discontinuity in the normal component, which is the sum of the normal viscous stress $\displaystyle\sigma_{n}$ and the pressure, is balanced by the surface tension

[TABLE]

where $n$ is a unit normal vector to the interface and double square brackets denote a jump in the corresponding quantity across the interface. The Weber number is defined as $\displaystyle\mbox{We}=\rho_{2}DU_{m}^{2}/\sigma$. Expressions for the stress tensor components are given in Appendix B.

3 Base flow

The base flow is assumed to be steady, laminar, and fully-developed flow. The only non-zero component is in the pipe axis direction, and it varies only in the cross section \displaystyle{\mbox{\boldmathU^{(k)}}}=\left(0,0,U_{z}^{(k)}(\xi,\phi)\right). Then the base flow momentum equations are obtained from Eqs. 7a-c and read

[TABLE]

Note that differentiation of Eqs. 13b,c in the axial $z$ direction, i.e., $\displaystyle\frac{\partial}{\partial z}$(Eqs. 13b,c), yields $\displaystyle\frac{\partial}{\partial\phi}\left(\frac{\partial P^{(k)}}{\partial z}\right)=0$ and $\displaystyle\frac{\partial}{\partial\xi}\left(\frac{\partial P^{(k)}}{\partial z}\right)=0$, therefore the base-flow pressure gradient in the axial direction is the same for the two fluids $\displaystyle\left(\frac{\partial P_{1}}{\partial z}=\frac{\partial P_{2}}{\partial z}=\frac{\partial P}{\partial z}\right)$.

The base flow velocity is subject to no-slip conditions at the pipe wall

[TABLE]

and continuous across the fluid–fluid interface

[TABLE]

while the tangential shear stress of the base flow is continuous as well

[TABLE]

For a plane unperturbed interface, continuity of the normal stress results in continuity of the pressure across the interface

[TABLE]

Integration of the base-flow velocity profiles over the lower- and upper-layer cross-sectional areas yields a relative input of each fluid into the prescribed total volumetric flow rate, which can be written as a ratio of the corresponding superficial velocity to the mixture velocity

[TABLE]

4 Linear stability

In the following, we study linear stability of above plane-parallel flow with respect to infinitesimal perturbations. The velocity and pressure of the perturbed flow are

[TABLE]

while the perturbed interface is described as a surface $\displaystyle\phi=\pi+\eta\left(\xi,z,t\right)$, where $\eta$ denotes the deformation of the interface from the undisturbed plane state. The infinitesimal perturbations are

[TABLE]

where $\alpha$ is the wavenumber of the perturbation.

To derive the linearized form of the boundary conditions B.13 at the unperturbed interface $\phi=\pi$, we use the Taylor expansion of a scalar function $f$ of the radius vector $r$ around $r_{0}$ in bipolar cylindrical coordinates

[TABLE]

where \displaystyle{\mbox{\boldmathr}}=\bigl{(}\xi,\pi+\eta\left(\xi,z,t\right),z\bigr{)}, \displaystyle{\mbox{\boldmathr_{0}}}=\bigl{(}\xi,\pi,z\bigr{)}, and \displaystyle{\mbox{\boldmathr}}-{\mbox{\boldmathr_{0}}}=\eta\left(\xi,z,t\right)H_{\phi}{\mbox{\boldmathe_{\phi}}} for infinitesimally small $\eta$. The function $f$ denotes either the axial velocity $u_{z}^{(k)}$ or the pressure $p^{(k)}$ in the derivations below.

4.1 Linearized governing equations

Linearizing the continuity 6 and momentum 7 equations for the perturbed flow and taking into account the corresponding equations for the base flow (Eq. 13), we get equations formulated for the perturbation amplitude (Eq. 21)

[TABLE]

[TABLE]

4.2 Linearized boundary conditions

In this section, we list the full set of the linearized boundary conditions. The no-slip condition at the pipe wall is formulated for the perturbation amplitude (Eq. 21) as

[TABLE]

The interfacial boundary conditions that were formulated in sec. 2.4 are linearized using the Taylor expansion 22 around $\phi=\pi$. The kinematic boundary condition (Eq. 10) for the perturbed flow reads

[TABLE]

The continuity of velocity components yields

[TABLE]

The interfacial boundary conditions require also continuity of the tangential stresses in planes \bigl{(}\xi,\phi\bigr{)} and \bigl{(}\phi,z\bigr{)} and a jump in the normal shear stress (Eqs. B.13a-c, respectively), and their linearized versions read

[TABLE]

5 Numerical method

5.1 Computational grid

The stability problem formulated in the previous section is solved on a staggered grid using the finite–volume method. The computational grid is represented by a pipe cross section divided into four-sided cells, whose faces coincide with the coordinate lines in the bipolar coordinates, so that each cell is rectangular as well as the whole domain.

We use the staggered grid, so that the scalar value of pressure $p$, the base flow velocity $U_{z}^{(k)}$, and the axial velocity of the perturbation $u_{z}^{(k)}$ are calculated in the cell centers defined by integer indices \big{[}\xi_{i},\phi_{j}\big{]}, where $i=0,1,...,N_{\xi}$ and $j=0,1,...,N_{\phi}$. For the sake of conciseness, we omit the tilde marks placed above the perturbation amplitude. The two other components of the perturbation velocity, $u_{\xi}^{(k)}$ and $u_{\phi}^{(k)}$, are calculated on the cell faces \displaystyle\big{[}\xi_{i+1/2},\phi_{j}\big{]} and \displaystyle\big{[}\xi_{i},\phi_{j+1/2}\big{]}, respectively, where \displaystyle\xi_{i+1/2}=\big{(}\xi_{i}+\xi_{i+1}\big{)}/2 and \displaystyle\phi_{j+1/2}=\big{(}\phi_{j}+\phi_{j+1}\big{)}/2. In addition, all the variables are calculated at the interface, which is aligned with the boundary between two rows of cells at $\phi=\pi$. Since the holdup is not known a priori and $\phi_{0}$ depends on it (see sec. 2.2), the mesh distribution in the $\phi$-direction is changing in the process of the base flow solution until the holdup, corresponding to the specified input flow rates of the phases, is found. The number of cells in the $\phi$-direction in each sublayer (phase) is set to be proportional to its relative height ($N_{\text{lower}}/N_{\phi}=(\pi-\phi_{0})/\pi$).

For the purpose of the finite–volume discretization of the governing equations, described below, the arc lengths of each cell have to be introduced (Issa,, 1988)

[TABLE]

Note that some other arc lengths are also required for the discretization, and they are defined as, e.g., $\displaystyle\Delta\xi_{i+1/2,j}=\int_{\xi_{i}}^{\xi_{i+1}}H_{\xi}d\xi$ and $\displaystyle\Delta\phi_{i,j+1/2}=\int_{\phi_{j}}^{\phi_{j+1}}H_{\phi}d\phi$, while $\displaystyle\Delta\xi_{i\_i-1/2,j}=\int_{\xi_{i-1/2}}^{\xi_{i}}H_{\xi}d\xi$ and $\displaystyle\Delta\phi_{i,j\_j-1/2}=\int_{\phi_{j-1/2}}^{\phi_{j}}H_{\phi}d\phi$. Their values are calculated in the same manner as $\displaystyle\Delta\xi_{i,j}$ and $\displaystyle\Delta\phi_{i,j}$ above.

The uniform cell distribution in the $\xi$-direction leads to very strong clustering of the cells near the triple points in the physical space, since $\Delta\xi_{i,j}\to 0$ for $\xi\to\pm\infty$ (Eq. LABEL:Eq:_Delta_XI). On the other hand, $\Delta\xi_{i,j}$ is relatively large near the pipe midplane resulting in relatively coarse grid. Therefore, we use hyperbolic tangent stretching

[TABLE]

where $\xi$ is originally uniformly distributed grid cells and $a_{\xi}$($=5$) is the stretching parameter that determines the degree of clustering, so that the cells are redistributed to be more dense near the midplane (around $\xi=0$). This can be seen in Fig. 2 for the cells along the interface (cf. minimal and maximal values of $\Delta\xi$ for red squares and blue dots). In the $\phi$-direction, the grid is stretched near the pipe walls and at both sides of the interface using the $sin$ function

[TABLE]

In this work we use $a_{\xi}=5$ and $a_{\phi}=0.12$.

5.2 Finite–volume base flow equations

The finite–volume form of the base flow equation 13a is obtained by integration of the differential equation over a computational cell with \displaystyle\xi\in\big{[}\xi_{i-1/2},\xi_{i+1/2}\big{]} and \displaystyle\phi\in\bigl{[}\phi_{j-1/2},\phi_{j+1/2}\bigr{]}. By multiplying the result of integration by $(\Delta\xi\Delta\phi)_{i,j}$ we get (c.f. Issa,, 1988)

[TABLE]

5.3 Finite–volume continuity equation

The finite–volume form of the continuity equation 23 is obtained by integration of the differential equation over a computational cell with \displaystyle\xi\in\big{[}\xi_{i-1/2},\xi_{i+1/2}\big{]} and \displaystyle\phi\in\bigl{[}\phi_{j-1/2},\phi_{j+1/2}\bigr{]}

[TABLE]

5.4 Finite–volume momentum equations

Integration of the momentum equation in $\xi$-direction 24a over \displaystyle\xi\in\big{[}\xi_{i},\xi_{i+1}\big{]} and \displaystyle\phi\in\bigl{[}\phi_{j-1/2},\phi_{j+1/2}\bigr{]} yields

[TABLE]

Integration of the momentum equation in the $\phi$-direction 24b over \displaystyle\xi\in\big{[}\xi_{i-1/2},\xi_{i+1/2}\big{]} and \displaystyle\phi\in\bigl{[}\phi_{j},\phi_{j+1}\bigr{]} yields

[TABLE]

Integration of the momentum equation in the $z$-direction 24c over \displaystyle\xi\in\big{[}\xi_{i-1/2},\xi_{i+1/2}\big{]} and \displaystyle\phi\in\bigl{[}\phi_{j-1/2},\phi_{j+1/2}\bigr{]} yields

[TABLE]

5.5 Finite–volume interfacial boundary conditions

The finite–volume discretization of the non-trivial (linear) boundary conditions at the interface (Eqs. 26, 27c, 28a-c), where all the variables are calculated at the boundary between two rows of cells at $\phi_{j}=\pi$. An infinitesimal interface deformation in the $y$-direction (in phisycal space) is defined through the Lamé coefficient

[TABLE]

Then the finite–volume form of the kinematic boundary condition (Eq. 26) is

[TABLE]

and the jump in the axial velocity of the perturbation at the interface (Eq. 26c) can be written as

[TABLE]

The finite–volume discretization of the linearized boundary conditions on the tangential stresses in planes \bigl{(}\xi,\phi\bigr{)} and \bigl{(}\phi,z\bigr{)} (Eqs. 28a-c)) yields

[TABLE]

Due to the surface tension, the pressure field is discontinuous across the interface. Therefore, at the interface the pressure in each phase has to be extrapolated from the neighboring cells

[TABLE]

and the finite–volume discretization of 28c reads

[TABLE]

5.6 Linear stability analysis

The discretization of the linearized governing equations (Eqs. 33 - 36) and boundary conditions(Eqs. 25 - 28c) reduces the linear stability problem to a generalized eigenvalue problem

[TABLE]

with the complex time increment $\lambda=\lambda_{R}+i\lambda_{I}$ being the eigenvalue and the perturbation \displaystyle{\mbox{\boldmathv}}=\big{(}u_{\xi},u_{\phi},u_{z},p,\Delta\eta\big{)}^{T} being the corresponding eigenvector. $B$ is a diagonal matrix with the elements corresponding to time derivatives of $\tilde{u}$ and $\eta$ equal to one. On the other hand, the diagonal elements corresponding to $p$ and the boundary conditions without time derivatives (all except for the kinematic one) are zeros, so that \det{\mbox{\boldmathB}}=0. The real part of the eigenvalue, $\lambda_{R}$ determines the growth rate of the perturbation, so that the flow is considered unstable when there exists an eigenvalue with a positive real part. Considering all possible wavenumbers, the eigenvalue with the largest $\lambda_{R}$ and the corresponding eigenvector are referred as leading or the most unstable one.

The generalized eigenvalue problem 43 is solved by the Arnoldi iteration in the shift-and-inverse mode

[TABLE]

where \displaystyle\vartheta=1/\bigl{(}\lambda-\lambda_{0}\bigr{)} and $\lambda_{0}$ is a complex shift. Following the approach of (Gelfgat,, 2007) and in similar manner as in (Gelfgat and Brauner,, 2020), we use the ARPACK package of (Lechouq et al.,, 1998) in FORTRAN for the Arnoldi iteration, computing the LU decomposition of the complex matrix \displaystyle({\mbox{\boldmathJ}}-\lambda_{0}\cdot{\mbox{\boldmathB}})^{-1}, so that the calculation of the next Krylov vector for the Arnoldi method is reduced to one backward and one forward substitution.

6 Immersed boundary method

An alternative numerical approach to solve the stability problem for two-phase stratified flow in a circular pipe is to use Cartesian formulation same as in Gelfgat and Brauner, (2020) and to treat the pipe wall by the immersed boundary method (IBM). In this approach, the computational domain is chosen as a square with the length of the side equals to the pipe diameter. Momentum forcing and mass source are applied outside of the flow domain and the velocity is interpolated across the pipe wall to satisfy no-slip boundary conditions. For this purpose, we apply interpolations of Kim et al., (2001) for the conservation of mass in the continuity equation and those of Wu, (2019) for velocities in the momentum equations. In the following section, the results obtained by IBM are cross verified with those in the bipolar coordinates.

7 Results and discussion

To the best of our knowledge, a rigorous linear stability analysis of two-phase stratified flow in a circular pipe has never been performed before. Since there are no independent results for comparison, we use two independent numerical approaches to verify our results. In addition, we can verify the numerical solution for the base flow by comparison with the exact analytical solution reported in Goldstein et al., (2015). Since the considered problem involves too many governing parameters, we focus on two representative examples of two-phase flows to demonstrate the proposed numerical approach. The physical properties and pipe diameters of the considered flows are summarized in Table 1.

7.1 Liquid–liquid pipe flow

As a benchmark liquid-liquid system, an oil-water flow in a circular pipe of diameter $D=0.02$m is chosen (see the first line in table 1). The light upper phase is oil with density and viscosity slightly different from water, $\rho_{12}=1.198$ and $\mu_{12}=0.103$. Once the physical properties and pipe diameter are defined, the holdup and base flow velocity profile are uniquely determined by the flow rate ratio (Goldstein et al.,, 2015). In this work, we evaluate the base flow by solving numerically equation 32 in the same manner as in Gelfgat and Brauner, (2020). In table 2, we demonstrate convergence of $U_{1S}/U_{2S}$ for two different holdups $h$ using the proposed numerical methods in the bipolar (denoted as bipolar FVM in the table) and Cartesian (denoted as IBM) coordinates. It can be seen that $U_{1S}/U_{2S}$ converges to the third decimal place (0.379) for the computational grids with 200 and more cells in each direction, for both numerical methods ($N_{\xi}(N_{x})=200,N_{\phi}(N_{y})=200$). The analytical solution of Goldstein et al., (2015) is also reported in table 2, and it can be seen that deviation of the numerical solution from the analytical one is much less than one percent.

The base flow velocity $U_{z}$ profile in the pipe cross-section is shown in Fig. 3. It is important to note that for $h=0.202$, although the lower (water) phase occupies only about a fifth of the pipe cross-sectional area as defined by its holdup, the height of the interface equals to about a fourth of the pipe diameter. Symmetry of the problem with respect to the midplane, $x=0$, prescribes symmetry of the base flow velocity that has only axial non-zero component. The base flow velocity at the liquid-liquid interface (horizontal dashed black line in Fig. 3a) is a function of the horizontal coordinate $x$ and reaches the maximum at the midplane, $x=0$, from which it decays to zero at the triple points (Fig. 3b). The overall maximum of the base flow velocity is located at the midplane below the interface in the bulk of water phase (Fig. 3a,c). These graphs also show that for the chosen grids there is no visible difference between the base flow velocity profiles for the two numerical methods.

Once the base flow solution is obtained, the critical operational conditions are found by changing the superficial velocities of both phases, while keeping their ratio constant, so that the holdup remains constant. The critical superficial velocities are defined as those for which the perturbations of all wavenumbers are damped ($\displaystyle\lambda_{R}<0$), except for a critical one $\displaystyle\alpha_{cr}$, for which the perturbation is neutrally stable, i.e., $\displaystyle Re[\lambda(\alpha_{cr})]=0$. For a holdup $h=0.202$, the critical wavenumber is found to be $\alpha_{cr}=5.6$. In table 3, we report the critical superficial velocities and their grid convergence. The critical conditions are also characterized by the imaginary part of the critical eigenvalue (see its convergence in table 3), which provides the oscillation frequency, $\displaystyle|\lambda_{I}|$, and the propagation speed of the critical perturbation, $c_{cr}=-\lambda_{I}/\alpha_{cr}$.

The results of the stability analysis are obtained on the computational grid with 200 cells in each direction, i.e., $\displaystyle N_{\xi}=N_{\phi}=200$, using the stretching mentioned in sec. 5.1. Such grid ensures convergence within three decimal places (table  3). Additionally, Fig. 4 illustrates convergence of the most unstable perturbation of the interface displacement (Fig. 4a) and of a jump in the pressure field across the interface due to the surface tension (Fig. 4b). An additional advantage of the bipolar coordinates compared to the Cartesian ones is a possibility to calculate flow fields very close to the triple points and to prescribe exactly the boundary conditions there, i.e., at $\xi\to\pm\infty$. The cells adjacent to the triple points have their center at $\pm\xi_{\text{max}}$. The base flow and critical eigenvalues become independent of $\xi_{\text{max}}$, when it is greater than $2\pi$ as demonstrated in table 4. Same convergence behavior in bipolar coordinates is found for the critical perturbation profiles, as shown for the interface displacement (a) and for the pressure jump at the interface (b) in Fig. 5. Note that the black dots ($\xi_{\text{max}}=3\pi$) in the bipolar coordinates in Fig. 5 correspond to the blue dots in Cartesian coordinates in Fig. 4. Contrarily to the base flow velocity reaching the maximum at the midplane, $x=0$, the interface displacement reaches its maximal value at $\xi\approx\pm 2$ that corresponds to $x\approx\pm 0.7$, while the maximum jump in pressure is attained closer to the triple point, at $\xi\pm 3.2$ ($x\approx\pm 0.8$). In this work, we use $\xi_{\text{max}}=3\pi$ (distance to the triple point $e^{-3\pi}\approx 8\times 10^{-5}$), so that the smallest grid size is $\approx 2\times 10^{-5}$. For comparison, the smallest grid size near the wall is approximately $0.005$ in the Cartesian grid (IBM) with $N_{x}=N_{y}=200$ and $sin$-stretching next to the walls (similar to Eq. 31). The finite distance from the triple point ($\approx e^{-\xi_{\text{max}}}$) necessarily introduces an error that is reflected in the values of $\lambda_{R}$ (table 4). At the instability threshold, these values are expected to be numerical zeroes (as for $\xi_{\text{max}}=3\pi$), so that the other values imply numerical errors. Far from the instability threshold, $\lambda_{R}$ converges similarly to the holdup and $\lambda_{I}$.

The velocity perturbation is three-dimensional, while its amplitude depends only on the two cross-sectional coordinates. Its amplitude is defined up to multiplication by a complex constant (see Eq. 21). Therefore, for the purpose of presentation of the critical perturbation, we scale the amplitudes by the maximal absolute value of the axial perturbation velocity max\big{(}|u_{z}|\big{)}, as is shown in Fig. 6. The amplitude of the axial velocity component is the largest (note that the limits of color bars, $\approx 0.45$ and $\approx 0.36$ in Fig. 6a,b, respectively, are less than one). Its two maxima are in the heavy phase between the midplane (vertical dash-dot black line) and the pipe wall Fig. 6a. The horizontal component, $u_{x}$, reaches the maximum at the midplane, while two secondary maxima (orange contours) develop slightly below the interface at $x\approx\pm 0.007$. The maxima of the vertical component of velocity, $u_{y}$, which has the smallest amplitude among the velocity components, is located at the interface at $x\approx\pm 0.008$.

Isolines of the most unstable perturbation obtained by the IBM with the grid of two hundred cells in each direction are indistinguishable from those shown in Fig. 6. The similarity of the perturbation profiles at the interface is shown in Fig. 7 for the horizontal and vertical components of velocity (a), the interface displacement (b), and the pressure difference across the interface (c). The discretization of the interface boundary conditions by the two numerical methods is noticeably different, however it does not affect the solutions behavior inside the pipe, except for the pressure jump in the vicinity of the triple points (Fig. 7c).The latter is nevertheless more accurately represented by the solution in the bipolar coordinates.

7.2 Gas–liquid pipe flow

In this section we consider an air-water flow in a circular pipe of diameter $D=0.014$m (line 2 of table 1). The air-water flow is characterized by large density and viscosity ratios, $\rho_{12}=1000$ and $\mu_{12}=55$, respectively. Compared to the oil-water flow considered above, a much lower water-to-air flow rate ratio is required for the same holdup of the heavy phase. For example, for a holdup $0.6$, the ratio of the superficial velocities of water to air $U_{1S}/U_{2S}=0.15$. The dimensionless base flow velocity in the pipe cross-section is shown in Fig. 8. As before, the velocity at the air-water interface (horizontal dashed black line in Fig. 8a) depends on the horizontal coordinate $x$ and is symmetrical around the midplane ($x=0$, Fig. 8b), where it reaches its maximum in the bulk of the light (air) phase (Fig. 8a,c). Due to the large water-to-air viscosity ratio, the water layer is perceived almost as a solid wall by the upper phase (air), except for the nonzero interfacial velocity of the water. There, the water reaches its maximal velocity, which is much lower than the maximum air velocity (compare $U_{z}(y=0)\approx 0.5$ with $\max(U_{z})\approx 4.3$). The air velocity profile is similar to that of the single-phase air flowing through the part of the pipe cross section it occupies. In horizontal air-water flows, a similar picture can be observed for the base flow for other water holdups (not shown).

While the short-wave perturbation ($\alpha_{cr}=5.6$) was found to be critical for the oil-water flow discussed in the previous section, for the considered air-water pipe flow, the long-wave perturbation ($\alpha_{cr}\to 0$) is critical for a holdup of $0.6$. To find small, but nonzero, value of the wavenumber $\alpha$ (i.e., large but finite wavelength), we compute the critical superficial velocities and the propagation speed $c=-Im(\lambda)/\alpha$ of the critical perturbation for several small $\alpha$ as summarized in table 5. In the following, the long-wave instability can be obtained starting from $\alpha=0.001$, from which the convergence up to the third decimal place for decreasing $\alpha$ is reached in the zero-wavenumber limit. The convergence of results for computational grids with different number of cells for FVM in the bipolar coordinates and with the IBM in the Cartesian coordinates are summarized in table 6. The grid convergence of the air-water flow is similar to that of the oil-water flow, and $N_{\xi}=N_{\phi}=200$ is sufficient for convergence of the critical superficial velocities and perturbation frequency. The convergence of perturbation profiles is shown in Fig. 9 for the interface displacement (a) and the horizontal (b) and vertical (c) velocities at the interface. Even though the components of velocity in the cross section have small absolute values compared to the axial velocity (as discussed below), the grid convergence can be clearly seen in the graphs.

One of the features of long-wave perturbation is that its axial component is dominating and is larger by several orders of magnitude than the two other components (Fig. 10, cf. the color bar limits). Moreover, for the long-wave perturbations, the amplitudes of the horizontal and vertical components scale down with $\alpha$ for long waves, so that $|u_{x}|/\alpha$ and $|u_{y}|/\alpha$ are independent on $\alpha\ll 1$. The maximum of the axial velocity (red contours) is at the center of the interface, and the flow on the air side is perturbed stronger, as depicted in Fig. 10c. Another local maximum in the air phase is located at the midplane above the base flow maximum (compare Fig. 10c with Fig. 8). Since the maximal value of the velocity field of the critical perturbation is at the interface, one can expect the onset of instability triggered by growth of the interface deformation. This instability can be classified as interfacial instability as opposed to the shear instability that was examined above for the oil-water flow. Further nonlinear analysis is required to find out whether such instability leads to transition to stratified wavy or other flow pattern. More details on the perturbation flow field at the interface can be seen in Fig. 11a (horizontal and vertical velocities), while the interface displacement is shown in Fig. 11b. These graphs show that the flow fields obtained by the two numerical methods are found to be similar.

8 Conclusions

Linear stability analysis of two-phase stratified flow in a horizontal circular pipe was studied numerically for the first time. The analysis took into account all infinitesimal three-dimensional disturbances, including disturbances of the interface shape. To avoid difficulties connected with cylindrical axis and definition of the interface position in commonly used Cartesian or cylindrical coordinates, the problem was treated in bipolar cylindrical coordinates. In these coordinates, both the pipe wall and the plane interface coincide with the coordinate surfaces. This allows us to apply the no-slip boundary conditions explicitly and impose the interfacial boundary conditions relatively easy. Additionally, it enables investigation of the local behavior of the flow fields and shear stresses in the vicinity of the triple point.

The results obtained in the bipolar coordinates were verified by an independent numerical solution based on the problem formulation in the Cartesian coordinates and the immersed boundary method. In this approach, the pipe wall is approximated by the immersed boundary technique, which slightly reduces the accuracy, while the plane unperturbed fluid-fluid interface remains a coordinate surface. The results computed by the two approaches are in good agreement, which verifies their correctness.

Two representative examples of gas-liquid and liquid-liquid flows with particular holdups were considered to demonstrate applicability of the proposed numerical technique. These cases were chosen such that we could verify the numerical solution for the base flow with the analytical solution derived in Goldstein et al., (2015), which was obtained in the same bipolar coordinate system.

In this work, we report results on the critical operational conditions and perturbations and show that our numerical approach accounting for all-wavelength perturbations, ranging from short- to very long-waves, allows computation of the shear (as in case of the air-water flow considered) and interfacial instabilities, as for the oil-water flow. Further comprehensive parametric studies, including computation of the stability boundaries for the entire range of flow holdups, are to be addressed in the future. An additional issue that is out of the scope of the present work is the case of nonplane unperturbed interface (e.g., constant-curvature interface in capillary-dominated flows, see Gorelik and Brauner,, 1999). Since the constant-curvature interface aligns with coordinate line of bipolar coordinate system, the proposed numerical method can be easily adjusted to such flows. In light of a recent study of Goldstein et al., (2021) on the local flow behavior near the triple point, future works should also address the effect of boundary conditions at the triple points on stability of two-phase stratified flow.

Acknowledgments

This research was supported by Israel Science Foundation (ISF) grant No 415/18.

Appendix A Governing equations in orthogonal coordinates

An orthogonal coordinate system $(\xi,\phi,z)$ is defined such that it is curvilinear in the cross-sectional coordinates $\xi$ and $\phi$, while $z$ is an axis coordinate identical to that in Cartesian coordinate system. In these coordinates, the velocity vector can be written through its coordinate components as \displaystyle{\mbox{\boldmathu}}^{(k)}=u_{\xi}^{(k)}{\mbox{\boldmathe_{\xi}}}+u_{\phi}^{(k)}{\mbox{\boldmathe_{\phi}}}+u_{z}^{(k)}{\mbox{\boldmathe_{z}}}, where \displaystyle\left\{{\mbox{\boldmathe_{\xi}}},{\mbox{\boldmathe_{\phi}}},{\mbox{\boldmathe_{z}}}\right\} form the orthonormal basis. The scale factors, also known as Lamé coefficients, are:

[TABLE]

and $H_{z}=1$. The vector equations 1 and 2 can be rewritten using the expressions of vector differential operators in orthogonal coordinates (Kochin et al.,, 1964):

[TABLE]

[TABLE]

Appendix B Interfacial boundary conditions

The fluid–fluid interface is defined as a surface

[TABLE]

In the present work, we considered a case of the plane unperturbed interface, i.e., $\displaystyle\phi^{*}=\pi$. The interface between the fluids is defined so that there is no mass transfer (i.e., fluid flow) across it. Therefore, the kinematic boundary condition at the interface, $\displaystyle\phi=\phi_{\eta}$, requires the material derivative of the surface with respect to time to be zero

[TABLE]

Continuity of the velocity at the interface reads

[TABLE]

A unit vector normal to the surface is defined as

[TABLE]

so the mean curvature of the interface is

[TABLE]

Linearization of the unit normal vector (Eq. B.4) yields

[TABLE]

the linearized version of the mean curvature of the interface then reads

[TABLE]

The linearized unit vectors tangential to the surface and located in the planes $(\xi,\phi)$ and $(\phi,z)$ are, respectively

[TABLE]

[TABLE]

The interfacial boundary conditions require continuity of the tangential components of the viscous stress tensor, \displaystyle{\mbox{\boldmath\mathcal{T}}}=\frac{\mu_{k}\mu_{12}}{\mu_{1}}\left(\nabla{\mbox{\boldmathu}}+\left(\nabla{\mbox{\boldmathu}}\right)^{T}\right), while the discontinuity in the normal component is balanced by the surface tension

[TABLE]

The tangential and normal viscous stresses are

[TABLE]

The components of the viscous stress tensor in bipolar coordinates are

[TABLE]

To derive boundary conditions for the viscous stresses, we substitute the components of the stress tensor ,B.12 and the normal and tangential unit vectors into Eq. B.11 to obtain the final expressions:

[TABLE]