A Unifying Numerical Framework for the “Small-Slope” Based Core-Annular Flow Instability Models

Abstract

The snap-off is an instability phenomenon that takes place during the immiscible two-phase flow in porous media due to competing forces acting on the fluid phases and at the interface between them. Different theoretical approaches have been proposed for the development of mathematical models that describe the dynamics of a fluid/fluid interface in order to analyze the snap-off mechanism. The models studied here are based on the “small-slope” approach and were derived from the mass conservation and other governing equations of two-phase flow at pore scale in circular capillaries for pure and complex interfaces. The models consist of evolution equations; highly nonlinear partial differential equations of fourth order in space and first order in time. Although the structure of the models for each type of interface is similar, different numerical techniques have been employed to solve them. Here, we propose a unifying numerical framework to solve the group of such models. Such a framework is based on the Fourier pseudo-spectral differentiation method which uses the Fast Fourier Transform (FFT) and the inverse FFT (IFFT) algorithms. We compared the solutions obtained with this method to the results reported in the literature in order to validate our framework. In general, acceptable agreements were obtained in the dynamics of the snap-off.

1. Introduction

When a fluid displaces another in a porous media, different displacement mechanisms may occur depending on the local flow conditions. The displacement may occur in a piston-like regime, at which one phase is fully pushed by the other, or the displacing liquid may leave a thin film at the wall forming an interface between the two fluids. If the capillary forces are strong enough, the thin film attached to the wall may form a liquid collar and could break up the displacing phase into separate ganglia due to surface tension. This instability phenomenon is known as snap-off (Figure 1). The occurrence of snap-off is determined by competing forces acting on the fluid phases and at the interface between them: the capillary pressure difference that drives the wetting phase towards the constriction, the shear stress at the interface that resists the capillary driven flow and the squeezing process of the non-wetting phase [1].

The snap-off phenomenon is important to many industrial and technological areas such as biofluid dynamics, CO₂ geological sequestration, oil recovery, microfluidic devices and chemical engineering. This phenomenon has been studied both experimentally and theoretically for viscous (drops) and inviscid (bubbles) fluids. This paper is focused on the theoretical approach. According to Deng et al., [2] the theoretical snap-off studies typically have two goals: development of models which describe the evolution of the interface between the fluids until snap-off is achieved and the definition of a snap-off criterion. The interface evolution models are used to test snap-off criteria and quantify the snap-off time.

Different approaches have been proposed for the development of mathematical models that describe the dynamics of a fluid/fluid interface [3,4,5,6,7,8]. Here we focused on the models that describe the drop and bubble snap-off in circular capillaries by means of an evolution equation derived from the conservation of mass equation and the small-slope approach [9,10]. This approach consists in the improvement of the approximation of the Young–Laplace equation; that is, in the balance of the pressure forces that have their expression in the curvature of the interface. The small slope approach assumes small capillary and Reynolds numbers, which allow the flow in the capillary to be approximated as Poiseuillean flow. Furthermore, it is assumed that the radius of the pore wall varies to a small extent relative to the axial distance, so lubrication approximation applies [11]. The models derived under this approach have been solved with diverse numerical techniques.

Hammond [12] performs a nonlinear analysis based on lubrication theory in order to examine the core-annular flow instability in straight and circular capillaries. The analysis results in a nonlinear thin-film evolution equation. To solve this model, he used the Method of Lines (MOL), which consists of replacing the spatial derivatives with finite differences that generate a system of coupled ordinary differential equations (ODE) for the time evolution of the gas–liquid interface; the ODE system is solved with the Gear’s method.

Gauglitz and Radke [9,10] extended the Hammond’s evolution equation by making higher order corrections in the thin-film analysis; this resulted in the small-slope approximation. The small-slope equations have a more complex mathematical form compared to the thin-film ones; this is due to a better approximation of the Young–Laplace equation. To analyze the film breakup in both straight and constricted cylindrical capillaries, they implement the Galerkin Finite Element (GFE) method to solve their model. They used Hermite cubic basis functions, and the time derivatives are discretized with the Crank–Nicholson procedure. The resulting system of algebraic equations is solved by iterations with the Newton–Raphson method.

Zhang et al. [13] solved the Gauglitz and Radke (GR) model for straight capillaries by using the MATLAB “pdepe” solver. However, they had to do a reduction of order transformation in order to adapt the problem to the algorithm requirements. The “pdepe” algorithm solves the initial boundary value problems for parabolic-elliptic partial differential equations (PDE) in 1-D using the MOL, where the spatial variable is discretized according with Skeel and Berzins [14] with Petrov–Galerkin’s approach. This results in a system of ordinary differential equations that are solved using a variable-step, variable-order algorithm based on the numerical differentiation formulas of orders 1 to 5 [15].

Recently, Hoyer et al. [16] extended the GR model including elastic interfacial stress to analyze the behavior of the breakup process in complex interfaces. The interfacial elasticity makes the surface stress anisotropic, which has a direct effect on the pressure along the liquid film. For inelastic cases, the Hoyer and coworkers (HAC) model becomes a more accurate formulation since the GR model linearizes the logarithmic terms and discards the film thickness terms of an order greater than three (thin-film assumption). Such discrimination can generate important differences, for example, in snap-off time (close to a factor of 2) when the models are tested at high capillary numbers. To solve their model, suitable for the case of dispersed drop viscosity much lower than the continuous phase, they used an explicit Euler first order finite difference method, also known as Forward Euler (FE).

Additionally, Beresnev and Deng [17] developed an evolution equation that describes the temporal evolution in a fluid–fluid interface for arbitrary viscosities in the presence of base flow. The Beresnev and Deng (BD) model is solved with the MOL framework implemented in Mathematica [18]. Furthermore, they validate their model with Computational Fluid Dynamics (CFD) experiments in two commercial models: FLUENT and COMSOL. To discretize the governing equations of flow, FLUENT uses the finite volume element method, while COMSOL use finite elements; for tracking the interface position the first uses the volume-of-fluid method and the second one the level-set method. Note that the computing time in the case of BD model is in the order of minutes, while CFD experiments took at least 1 day per scenario.

It is from this general review, that we can appreciate that there is great diversity in the schemes traditionally used to solve the models of the different interfaces reported (gas–liquid elastic, gas–liquid and liquid–liquid). This situation may be considered as restrictive when it is necessary to assess the effects of different fluids in a two-phase flow system under instability conditions. We, therefore, propose a unifying pseudo-spectral differentiation framework in which the most important reported models are solved. To validate our framework, we compare the solutions obtained using this method with those reported by our reference authors utilizing other numeric techniques.

This article is divided into four parts. Section 2, called Materials and Methods, first presents the geometry of the problem. Later, some considerations to derive the evolution equations from the governing equations of flow at pore-scale are briefly presented. Since the formulation of each model has its own particularities that depend on the fluid/fluid system, the main assumptions considered to represent the rheological nature of each type of interface are also mentioned. Then the models are presented in their dimensionless form, as well as the corresponding boundary and initial conditions. The last part of this section is dedicated to present Fourier pseudo-spectral differentiation method. Following that, Section 3 includes the results of the evaluation of the models using the framework proposed here, as well as a comparison with those reported in the literature. Section 4 contains the conclusions of our research.

2. Materials and Methods

2.1. Geometry of the Problem

The porous structure of a reservoir presents complex geometric shapes; among these, the presence of constrictions in the capillary channels called pore throats is found. Constricted pore geometry is typical of natural core-annular flows and can determinate the fluid breakup occurrence. During core-annular flow, the pores have a preference to be wetted by one of the phases, either water or oil, whereby one layer of the wetting phase will be present between the pore-wall and the non-wetting phase. This fluids configuration and geometry are schematized in Figure 2.

This basic geometry of the porous medium, consisting of a circular capillary with one or more constrictions, has been used to study the mechanism of drop breakup during the immiscible two-phase flow at pore-scale. Although any smooth function can be used to represent this particular geometry, a sinusoidal function is typically used [17]:

$$\mathsf{\lambda }\left(x\right)=\frac{1}{2}{R}_c\left[\left(1+a\right)-\left(1-a\right)\cos \left(\frac{\mathsf{\alpha }\mathsf{\pi }x}{\left(1-a\right)R_c}\right)\right],$$

(1)

where $\mathsf{\lambda }$ is the radial coordinate of the capillary wall, $x$ refers to the axial coordinate, $a$ is the ratio between the throat ($R_t$) and capillary ($R_c$) radii and $\mathsf{\alpha }$ is a geometric parameter which represents the ratio between the capillary radius ($R_c$) and half of the capillary length $(\ell )$, $\mathsf{\alpha }=R_c/(\ell /2)$.

2.2. Governing Equations

The governing equations for a two-phase flow at pore scale are the simplified Navier–Stokes equations (momentum equation) using the lubrication approximation for incompressible and laminar flow in cylindrical coordinates, the conservation of mass, continuity and the Young–Laplace law derived under small-slope approach [17]:

$$-\frac{\partial P_{1,2}}{\partial x}+{\mathsf{\mu }}_{1,2}\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial u_{1,2}}{\partial r}\right)=0,$$

(2)

$$\left(\frac{\partial Q_1^{}}{\partial x^{}}\right)dx=-2\mathsf{\pi }{R}_c{}^2\mathsf{\kappa }\left(\frac{\partial \mathsf{\kappa }}{\partial t}\right)dx,$$

(3)

$$Q_1^{}+Q_2^{}=Q,$$

(4)

$$P_c=P_1-P_2=\sigma \left(\frac{1}{\mathsf{\kappa }}-\frac{{\partial }^2\mathsf{\kappa }}{\partial x^2}\right),$$

(5)

where $P$ is the pressure, subscripts 1 and 2 denote, respectively, the non-wetting and wetting phase, $\mathsf{\mu }$ is the dynamic viscosity, $r$ the radial coordinate, $u$ is the axial velocity, $t$ the temporal variable, $Q$ the total volumetric flow rate in the capillary, $P_c$ is the capillary pressure, $\mathsf{\kappa }$ the radius of interface and $\sigma$ the surface tension.

The drop breakup problem needs to determine the interface position change over time $\partial \mathsf{\kappa }\left(x,t\right)/\partial t$, thus the equations that govern the two-phase flow at pore scale are combined to establish an evolution equation that describes the interface dynamics. The evolution equation is derived from the mathematical expression of the conservation of mass (Equation (3)), which relates the change rate of the volume of some of the phases with the volumetric flow differential, which can be written in the form:

$$\frac{\partial \mathsf{\kappa }}{\partial t}=-\frac{1}{2\mathsf{\pi }{R}_c{}^2\mathsf{\kappa }}\left(\frac{\partial Q_1^{}}{\partial x^{}}\right),$$

(6)

In order to close the evolution equation, the volume flux of the non-wetting phase $Q_1$ is required to be explicitly expressed through the interface position $\mathsf{\kappa }$. The evolution equation derived from this process will depend on the nature of the fluid phases, though its own particularities are further detailed in Section 2.3.1.

2.3. Drop and Bubble Snap-Off Models

In this section, we present the reported model’s equations that describe the snap-off dynamics in circular capillaries for elastic gas–liquid, gas–liquid and liquid–liquid interfaces. Even if original works can be consulted for further details of mathematical procedure for models derivation [9,10,16,17], it is important to specify the physical assumptions taken for such derivation.

2.3.1. General Considerations for the Derivation of Models

In general terms, the evolution equations that describe the temporal evolution of the interface are basically obtained by solving the simplified incompressible Navier–Stokes equations (momentum equation) for cylindrical core-annular geometry (Equation (2)); with which the velocity profiles of the wetting and non-wetting phases are obtained. The velocity profiles are integrated to obtain the volume flow rates. The resultant formulations are combined with the continuity and Young–Laplace equations (Equations (4) and (5)) to obtain an explicit expression of the volume flux through the interface position $\mathsf{\kappa }$, which is finally substituted in the conservation equation (Equation (6)).

Depending on the core phase nature, two approaches can be distinguished: (1) those that consider an inviscid core fluid and constant zero pressure (GR and HAC models) (2) those in which a viscous core liquid is studied (BD model). For both cases, non-slip (solid wall) at the surface of the capillary ($r=\mathsf{\lambda }\hspace{0.17em}\to \hspace{0.17em}{u}_2=0$) and balance of forces at the interface $\left(r=\mathsf{\kappa }\hspace{0.17em}\to \hspace{0.17em}{\mathsf{\mu }}_{1}d{u}_1/dr={\mathsf{\mu }}_{2}d{u}_2/dr\right)$ are postulated as boundary conditions. Additionally, this model assumes that the interface is infinitesimally thin and in a state of physic-chemical and mechanical equilibrium; therefore, no transport phenomena occurs.

In the first approach, since the viscosity of the non-wetting phase is null (${\mathsf{\mu }}_1=0$), shear stress of the core is neglected ${\mathsf{\mu }}_{1}d{u}_1/dr=0$; consequently $r=\mathsf{\kappa }\hspace{0.17em}\to \hspace{0.17em}{\mathsf{\mu }}_{2}d{u}_2/dr=0$. In Hoyer and coworkers formulation [16] for elastic interfaces, surfaces stress ($T=\sigma +K\xi$) in axial ($xx$) and transversal ($\mathsf{\theta }\mathsf{\theta }$) directions are included as function of the surface tension $\sigma$, the dilational interfacial elasticity $K$ and the corresponding surface strains ${\xi }_{xx}$ and ${\xi }_{\mathsf{\theta }\mathsf{\theta }}$ (defined by geometric relationships). Note that for inelastic cases $\left(K=0\right)$, the surface stress becomes $T_{xx}=T_{\mathsf{\theta }\mathsf{\theta }}=\sigma$. At the interface the liquid film shear stress is balanced by the net surfaces stress in axial direction, then: $r=\mathsf{\kappa }\hspace{0.17em}\to \hspace{0.17em}{\mathsf{\mu }}_{2}d{u}_1/dr=d{T}_{xx}/dx$. Furthermore, a constant zero pressure is assumed for inviscid core ($P_1=0$), thus explicit expressions based on the interface curvature for the capillary pressure can be used: for inelastic case the Young–Laplace equation becomes $P_c=-P_2=\sigma \left[1/\mathsf{\kappa }+{\partial }^2\mathsf{\kappa }/\partial x^2\right]$; while for an elastic interface the expression is $P_c=-P_2=\left[\left(\sigma +K{\xi }_{\mathsf{\theta }\mathsf{\theta }}\right)/\mathsf{\kappa }+\left(\sigma +K{\xi }_{xx}\right)\left({\partial }^2\mathsf{\kappa }/\partial x^2\right)\right]$.

For fluid arbitrary viscosities (${\mathsf{\mu }}_1\ne 0$), the velocity profile of the core is equal to that of the annular film at the interface ($r=\mathsf{\kappa }\hspace{0.17em}\to \hspace{0.17em}{u}_1=u_2$); additionally no variation of the radial velocity derivative on the axes is postulated $\left(r=0\hspace{0.17em}\to \hspace{0.17em}d{u}_1/dr=0\right)$. In this case the assumption of constant zero pressure in the core is no longer valid, and Laplace’s law would only provide a relative pressure difference between the two fluids, it is expressed as: $P_c=P_1-P_2=\sigma \left[1/\mathsf{\kappa }+{\partial }^2\mathsf{\kappa }/\partial x^2\right]$.

The nondimensionalizing is performed according with:

$$\begin{gathered} L=\ell /R_c; x^{*}=x/R_c; {\mathsf{\lambda }}^{*}=\mathsf{\lambda }/R_c; {\mathsf{\kappa }}^{*}=\mathsf{\kappa }/R_c; {\delta }^{*}=\delta /R_c ; \\ {P}^{*}=P/\left(\sigma /R_c\right); \\ {Q}^{*}=Q/\left(\sigma R_c{}^2/{\mathsf{\mu }}_2\right); \mathsf{\tau }=t/\left(3{\mathsf{\mu }}_2{R}_c/\sigma \right); K^{*}=K/\sigma ; u^{*}=u/\left(\sigma /{\mathsf{\mu }}_2\right). \end{gathered}$$

(7)

where * denotes dimensionless variables and $\delta$ is the film thickness. For the viscous core case, the dimensionless equations use the viscosity of the non-wetting phase ${\mathsf{\mu }}_1$ instead of the wetting phase ${\mathsf{\mu }}_2$, and the dimensionless time is $\mathsf{\tau }=t/\left({\mathsf{\mu }}_1{R}_c/\sigma \right).$

2.3.2. Inviscid Core Fluid Models

Hoyer and coworkers (HAC) model for an elastic behavior of a gas–liquid interface, encapsulated in $K*=0$ parameter, is denoted by the following formulation [16]:

$$\begin{array}{l}\frac{\partial {\mathsf{\kappa }}^{*}}{\partial {\mathsf{\tau }}^{*}}=-\frac{3}{16{\mathsf{\kappa }}^{*}}\left\{\left(\frac{1}{{\mathsf{\kappa }}^{*}{}^2}\frac{\partial {\mathsf{\kappa }}^{*}}{\partial x^{*}}+\frac{{\partial }^3{\mathsf{\kappa }}^{*}}{\partial x^{*}{}^3}\right)\left[-4{\mathsf{\lambda }}^{*}{}^2{\mathsf{\kappa }}^{*}{}^2+{\mathsf{\lambda }}^{*}{}^4+4\ln \left(\frac{{\mathsf{\lambda }}^{*}}{{\mathsf{\kappa }}^{*}}\right){\mathsf{\kappa }}^{*}{}^4+3{\mathsf{\kappa }}^{*}{}^4\right]\right\}+\\ \frac{\partial }{\partial x^{*}}\left\{K^{*}\frac{1}{\sqrt{1+{\frac{\partial {\mathsf{\kappa }}^{*}}{\partial x^{*}}}^2}}\left\{\begin{array}{l}4{\mathsf{\lambda }}^{*}{}^2{\mathsf{\kappa }}^{*}\left\{-1+\left[1-2\ln \left(\frac{{\mathsf{\kappa }}^{*}}{{\mathsf{\lambda }}^{*}}\right)\right]{\left(\frac{{\mathsf{\kappa }}^{*}}{{\mathsf{\lambda }}^{*}}\right)}^2\right\}\left(\frac{\partial {\mathsf{\kappa }}^{*}}{\partial x^{*}}\frac{{\partial }^2{\mathsf{\kappa }}^{*}}{\partial x^{*}{}^2}\right)-\\ {\mathsf{\lambda }}^{*}{}^4\left\{1-4{\left(\frac{{\mathsf{\kappa }}^{*}}{{\mathsf{\lambda }}^{*}}\right)}^2+\left[3-4\ln \left(\frac{{\mathsf{\kappa }}^{*}}{{\mathsf{\lambda }}^{*}}\right)\right]{\left(\frac{{\mathsf{\kappa }}^{*}}{{\mathsf{\lambda }}^{*}}\right)}^4\right\}\\ \left[\frac{\partial {\mathsf{\kappa }}^{*}}{\partial x^{*}}{\left(\frac{{\partial }^2{\mathsf{\kappa }}^{*}}{\partial x^{*}{}^2}\right)}^2+\left(1+{\frac{\partial {\mathsf{\kappa }}^{*}}{\partial x^{*}}}^2\right)\frac{{\partial }^3{\mathsf{\kappa }}^{*}}{\partial x^{*}{}^3}\right]\end{array}\right\}\right\}\end{array}$$

(8)

if $K*=0$, i.e., for inelastic interface, it then results in an extension of the GR model [16]:

$$\frac{\partial {\mathsf{\kappa }}^{*}}{\partial {\mathsf{\tau }}^{*}}=-\frac{3}{16{\mathsf{\kappa }}^{*}}\frac{\partial }{\partial x^{*}}\left\{\left(\frac{1}{{\mathsf{\kappa }}^{*}{}^2}\frac{\partial {\mathsf{\kappa }}^{*}}{\partial x^{*}}+\frac{{\partial }^3{\mathsf{\kappa }}^{*}}{\partial x^{*}{}^3}\right)\left[-4{\mathsf{\lambda }}^{*}{}^2{\mathsf{\kappa }}^{*}{}^2+{\mathsf{\lambda }}^{*}{}^4+4\ln \left(\frac{{\mathsf{\lambda }}^{*}}{{\mathsf{\kappa }}^{*}}\right){\mathsf{\kappa }}^{*}{}^4+3{\mathsf{\kappa }}^{*}{}^4\right]\right\},$$

(9)

The Gauglitz and Radke model (GR) for the gas-liquid interface is [10]:

$$\frac{\partial {\mathsf{\kappa }}^{*}}{\partial {\mathsf{\tau }}^{*}}=-\frac{1}{{\mathsf{\kappa }}^{*}}\frac{\partial }{\partial x^{*}}\left\{\left(\frac{1}{{\mathsf{\kappa }}^{*}{}^2}\frac{\partial {\mathsf{\kappa }}^{*}}{\partial x^{*}}+\frac{{\partial }^3{\mathsf{\kappa }}^{*}}{\partial x^{*}{}^3}\right)\left[{\mathsf{\lambda }}^{*}{\left({\mathsf{\lambda }}^{*}-{\mathsf{\kappa }}^{*}\right)}^3\right]\right\},$$

(10)

here the second factor on the derivative corresponds to an approximation of the equivalent term in the HAC model. For straight capillaries, i.e., ${\mathsf{\lambda }}^{*}\left(x^{*}\right)=1$, the GR model results in:

$$\frac{\partial {\mathsf{\kappa }}^{*}}{\partial {\mathsf{\tau }}^{*}}=-\frac{1}{{\mathsf{\kappa }}^{*}}\frac{\partial }{\partial x^{*}}\left[\left(\frac{1}{{\mathsf{\kappa }}^{*}{}^2}\frac{\partial {\mathsf{\kappa }}^{*}}{\partial x^{*}}+\frac{{\partial }^3{\mathsf{\kappa }}^{*}}{\partial x^{*}{}^3}\right){\left(1-{\mathsf{\kappa }}^{*}\right)}^3\right],$$

(11)

2.3.3. Viscous Core Fluid Models

The Beresnev and Deng model (BD) for liquid–liquid interfaces is [17]:

$$\begin{array}{l}\frac{\partial {\mathsf{\kappa }}^{*}}{\partial {\mathsf{\tau }}^{*}}=-\frac{1}{2\mathsf{\pi }}\frac{{\mathsf{\mu }}_2}{{\mathsf{\mu }}_1}\frac{Q^{*}}{{\mathsf{\kappa }}^{*}}\frac{\partial }{\partial x^{*}}\left[\frac{1+2\frac{{\mathsf{\mu }}_1}{{\mathsf{\mu }}_2}\left(\frac{{\mathsf{\lambda }}^{*}{}^2}{{\mathsf{\kappa }}^{*}{}^2}-1\right)}{\frac{{\mathsf{\lambda }}^{*}{}^4}{{\mathsf{\kappa }}^{*}{}^4}+\frac{{\mathsf{\mu }}_2}{{\mathsf{\mu }}_1}-1}\right]\\ +\frac{1}{4{\mathsf{\kappa }}^{*}}\left\{{\mathsf{\kappa }}^{*}{}^2\left(\frac{\partial {\mathsf{\kappa }}^{*}}{\partial x^{*}}+{\mathsf{\kappa }}^{*}{}^2\frac{{\partial }^3{\mathsf{\kappa }}^{*}}{\partial x^{*}{}^3}\right)\left\{\frac{\left(\frac{{\mathsf{\lambda }}^{*}{}^2}{{\mathsf{\kappa }}^{*}{}^2}-1\right)\left[\left(\frac{{\mathsf{\mu }}_1}{{\mathsf{\mu }}_2}-\frac{1}{4}\right)\frac{{\mathsf{\lambda }}^{*}{}^2}{{\mathsf{\kappa }}^{*}{}^2}+\frac{3}{4}-\frac{{\mathsf{\mu }}_1}{{\mathsf{\mu }}_2}\right]\hspace{1em}}{\frac{{\mathsf{\lambda }}^{*}{}^4}{{\mathsf{\kappa }}^{*}{}^4}+\frac{{\mathsf{\mu }}_2}{{\mathsf{\mu }}_1}-1}-\frac{{\mathsf{\mu }}_1}{{\mathsf{\mu }}_2}\ln \frac{{\mathsf{\lambda }}^{*}}{{\mathsf{\kappa }}^{*}}\right\}\right\}\end{array}.$$

(12)

2.4. Initial and Boundary Conditions

The models presented here consist of fourth order partial differential equation in space and first order in time. A uniform film thickness $\delta$ is assumed as an initial condition in all cases. For the HAC and GR models, fourth boundary conditions are required, while in BD model periodic boundary conditions are postulated.

2.4.1. Initial Conditions

The initial conditions are related to the residual film of liquid wetting the adjacent wall of the capillary that is formed when it is displaced by the non-wetting fluid. When the invading phase is gas (HAC and GR models), Bretherton’s theory [19] predicts the dimensionless thickness of the wetting film as a power law of the capillary number $\left(Ca\right):$

$${\delta }^{*}=1.337C{a}^{2/3},$$

(13)

Beresnev et al. [20] developed a hydrodynamic theory and performed microfluidic experiments for two liquids core-annular flow (BD model). They concluded that Bretherton’s equation does not apply to liquid–liquid displacement. Deng et al. [21] used these experimental measurements and fitted them to propose an empirical equation in order to predict the film thickness under these flow conditions:

$${\delta }^{*}=0.0412lo{g}_{10}\left(Ca\right)+0.1475,$$

(14)

2.4.2. Boundary Conditions

For HAC and GR models, the two first boundary conditions are related with the assumption that the bubble is centered at the neck of the throat, consequently the problem is considered symmetric. Thus the interface profile should be perpendicular to the plane $x^{*}=0$ and no flow at the pore neck is imposed; i.e., the first and third derivatives are equal to zero at $x^{*}=0$:

$${\frac{\partial {\mathsf{\kappa }}^{*}\left(x^{*},\mathsf{\tau }\right)}{\partial x^{*}}|}_{x^{*}=0}=0,$$

(15)

$${\frac{{\partial }^3{\mathsf{\kappa }}^{*}\left(x^{*},\mathsf{\tau }\right)}{\partial x^{*3}}|}_{x^{*}=0}=0,$$

(16)

The rest of the boundary conditions are imposed away of the constriction, at a distance $d$, where no movement in the interface is detected, thus:

$${{\mathsf{\kappa }}^{*}\left(x^{*},\mathsf{\tau }\right)|}_{x^{*}=d}={\mathsf{\lambda }}^{*}\left(d\right)-\delta ,$$

(17)

$${\frac{\partial {\mathsf{\kappa }}^{*}\left(x^{*},\mathsf{\tau }\right)}{\partial x^{*}}|}_{x^{*}=d}={\frac{\partial {\mathsf{\lambda }}^{*}\left(x^{*}\right)}{\partial x^{*}}|}_{x^{*}=d},$$

(18)

Periodic boundary conditions for the BD model are imposed, that is, the interface radius is equal at the boundaries for all simulation time:

$${\frac{\partial {\mathsf{\kappa }}^{*}\left(x^{*},\mathsf{\tau }\right)}{\partial x^{*}}|}_{x^{*}=-L/2}={\frac{\partial {\mathsf{\kappa }}^{*}\left(x^{*},\mathsf{\tau }\right)}{\partial x^{*}}|}_{x^{*}=L/2}.$$

(19)

It should be noted that with this boundary conditions gas-liquid models describe a symmetric interface evolution that reduced the computing to the middle of the domain, while the liquid–liquid formulation does not specify symmetry, this implies a calculation on the entire domain.

2.5. Numerical Procedures

The models presented in Section 2.3.2 and Section 2.3.3 consist in highly non-linear partial differential equations, fourth order in space and first order in time, whose analytical solution can be very complex if no simplifications are made. To solve such equations several numerical techniques have been proposed, which range from the simple Forward Euler method to the use of powerful solvers implemented in commercial software. We have thus proposed a unifying framework based on Fourier pseudo-spectral differentiation method to solve the group of exposed models.

Fourier-Based Pseudo-Spectral (FPS) Differentiation Method

The finite differences differentiation method approximates the derivatives of a function in a node of interest, where the accuracy order is related to the numbers of nodes involved in the differentiation formula: two nodes for a second-order accuracy $O(h^2)$, four nodes for a $O(h^4)$, etc. From a matricial point of view, the accuracy order determines the density of the differentiation matrix: to a second-order accuracy $O(h^2)$ approximation corresponds a tridiagonal differentiation matrix, to an $O(h^4)$ pentadiagonal matrix. The spectral methods have their basis in taking this process to the limit, i.e., it works with differentiation matrix of N-order of accuracy $O(h^N)$, where N is the total number of nodes, although this implies a dense matrix [22].

In the spectral collocation method (also known as pseudo-spectral), the unknown solution to the differential equation is expanded as a global interpolant, such as trigonometric or polynomial interpolant [23]. Since the geometry of our problem consists of a periodic domain (sinusoidal) and periodic boundary conditions are used, it results adequate for the use of a global trigonometric interpolant based on the compute of the Discrete Fourier Transform (DFT) and the inverse DFT (IDFT).

The DFT of a periodic function, in the domain $x^{*}=\left[-L/2\hspace{1em}L/2\right]$ divided on $N$ nodes spaced for a distance $h=2\mathsf{\pi }/N,$ period $T=2\mathsf{\pi }$, wave numbers $k$ entire multiples of $\mathsf{\pi }/h=N/2$, is stablished by:

$$\widehat{\mathsf{\kappa }}*=h{\displaystyle \sum _{j=1}^N{e}^{\mathrm{i}k{x}_j}}\mathsf{\kappa }{*}_j,$$

(20)

for all $k=-N/2+1,-N/2,\dots ,N/2$; where $\widehat{\mathsf{\kappa }}$ is the DFT, $i$ is the imaginary unit and $k$ is the variable in the Fourier space or wave number. While the inverse Fourier transform is defined as:

$${\mathsf{\kappa }}_j*=\frac{1}{2\mathsf{\pi }}{\displaystyle \sum _{k=-N/2+1}^{N/2}{e}^{\mathrm{i}k{x}_j}}{\widehat{\mathsf{\kappa }}}_k,$$

(21)

To compute the DFT and IDFT we used the Cooley and Tukey [24] algorithms, known as Fast Fourier Transform (FFT) and the inverse FFT (IFFT) implemented on MATLAB. Figure 3 shows a diagram with a flowchart that illustrates the steps followed to perform the pseudo-spectral differentiation method [25]. The procedure used to obtain the n-order spatial derivative of the dimensionless interface position at time $\mathsf{\tau }$, $\mathsf{\kappa }{*}^{\prime }=\partial \mathsf{\kappa }\left(x*,\mathsf{\tau }\right)/\partial \mathsf{\tau }$, on the domain $x^{*}=\left[-L/2\hspace{1em}L/2\right]$, begins with the compute of the DFT by the FFT ($\widehat{\mathsf{\kappa }}*$). Subsequently the n-order derivative of $\widehat{\mathsf{\kappa }}*$ is performed on the Fourier space by using the associated wave numbers $k$. Finally, the derivative $\mathsf{\kappa }{*}^{\prime }$ is obtained via the IFFT.

Once obtained the n-order spatial derivatives of the interface position ${\mathsf{\kappa }}^{*}$ that intervene on the evolution equations (Equations (8)–(12)), the first order temporal derivative is determined using a forward finite difference, in such a way that the equation becomes explicit and the solution of next time step can be estimated.

3. Results and Discussion

This section has the objective of validating our unifying numerical framework. For this purpose, we evaluated some representative cases of the models studied here in order to compare the results obtained with the Fourier-based pseudo-spectral (FPS) differentiation method to those reported in the literature with other numerical techniques.

To solve the HAC model for gas–liquid elastic interfaces (Equation (8)) the authors used the Forward Euler (FE) method [16]. Here, we focused on those results that show the effect of the interfacial elasticity [26], encapsulated in the dimensionless parameter $K^{*}$. The HAC model is evaluated in a geometry defined by $a=0.2$, $\mathsf{\alpha }=0.1$. The initial condition is defined by $Ca=0.1$. First, Figure 4 shows the evolution of the interface radius ${\mathsf{\kappa }}^{*}$ over time at the center of the capillary for cases where $K^{*}$ values are (a) less than and (b) greater than the unit. In first case ($K^{*}=0.4$) snap-off occurs, but not in second case ($K^{*}=1.1$). While in Figure 5, it is related to the interface profiles moments before the end of the simulation for different values of $K^{*}$.

The solutions reported for the GR model [10] (gas–liquid interfaces), represented by Equation (10), were obtained using the GFE method, which are compared to those obtained by this work with the FPS method. Figure 6 shows the interface evolution profiles at different dimensionless times, for the same geometry used on the previous case but with initial conditions defined by $Ca=8.94\times {10}^{-4}$. In the case of a straight capillary, the GR model is evaluated for a dimensionless initial film thickness $\mathsf{\epsilon }=0.12$ and a dimensionless wavelength of fastest growing disturbance of ${\mathsf{\lambda }}_{\max }^{*}=2^{3/2}\mathsf{\pi }$. In Figure 7 we compared the evolution of the interface profiles reported by Gauglitz and Radke [9] using the Galerkin FE method, to those obtained with the MATLAB “pdepe” solver by Zhang et al. [13] and to the results of this work applying the FPS method.

Regarding the viscous core fluid case, the BD model for liquid–liquid interfaces (Equation (12)) is first evaluated for a base case defined by $Ca={10}^{-3}$, $a=0.6$, ${\delta }^{*}=0.5R_t{}^{*}$, ${\mathsf{\mu }}_1/{\mathsf{\mu }}_2=10$ and different values of $\mathsf{\alpha }$. This base case is used by the BD model authors to validate the numerical solutions obtained with the MOL implemented on Mathematica [18] through the comparison of the solutions obtained by numerical simulations performed on CFD commercial codes. Here, we compare the results of Beresnev and Deng [17] with those obtained on this work with the FPS method (see Figure 8). The snap-off time for the cases that appear in Figure 8 are: (a) ${\mathsf{\tau }}_s=3.35$ (b) ${\mathsf{\tau }}_s=7.74$ (c) ${\mathsf{\tau }}_s=39.54$. Beresnev and Deng (Figure 3 in [14]) reported ${\mathsf{\tau }}_{sBD}=3.4,\hspace{0.17em}7.8,\hspace{0.17em}38.0$, respectively.

After these comparisons, the BD model is evaluated for two specific cases depicted in Figure 9. The first inspects a geometric break-up criterion (Equation (6) in [23]); the geometric slope overcomes the upper limit of such criterion ($\mathsf{\alpha }=0.729$) with the rest of parameters maintained. This means that, in theory, the snap-off is inhibited. In one of the profiles shown in Figure 9 it may be appreciated that the FPS reproduces this situation, i.e., snap-off does not occur and the wetting liquid film decreases its thickness. The second case in Figure 9 corresponds to a reversed viscosity ratio, that is the most viscous liquid is now the wetting phase, thus: ${\mathsf{\mu }}_1/{\mathsf{\mu }}_2=0.1$. The rest of the parameters are those of the base case, but with $Ca=3\times {10}^{-5}$ and $\mathsf{\alpha }=0.0695$; this last parameter corresponds with the geometry in case of Figure 8c. In fact, this case tries to demonstrate that even in most unfavorable conditions for drop breakup (displaced wetting phase more viscous than the non-wetting phase and a low flow regime), the breakup occurs if the geometrical criterion of rupture is met. Although, it should be noted that the effect of these “unfavorable” flow conditions is reflected on the snap-off time ${\mathsf{\tau }}_s$; since this is increased from ${\mathsf{\tau }}_s=39.54$ to ${\mathsf{\tau }}_s=600.16$. Beresnev and Deng (Figure 3c and Figure 7 in [14]) report for these cases ${\mathsf{\tau }}_{sBD}=38.0,\hspace{0.17em}585.0$, respectively.

As it may be observed on the series of Figure 4, Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9, there is an acceptable agreement between the results of the snap-off dynamics obtained in this work and those reported on the literature. The small deviations on snap-off time and form of the interface can be related with the arbitrary definition of breakup time. Regarding the differences in the interface evolution of inviscid and viscous core fluid, it should be noted that in the first case the forms of the interface are axisymmetric, in addition to the fact that a single breaking point is observed (Figure 5 and Figure 6). While in the second case the results show a more varied interface dynamics, both in non-symmetrical forms and in the presence, in some configurations, of more than one breaking point (Figure 8 and Figure 9). This can be explained by the high non-linearity of the BD model (Equation (12)) and the periodic boundary conditions imposed (Equation (19)).

It is important to remark that the small-slope based models assume $\mathsf{\alpha }\ll 1$, so the smaller $\mathsf{\alpha }$, the better the approximation. On Beresnev et al. (Table III in [23]), for example, it is clear that the order of the ratio between the theoretical and experimental snap-off times decreases for greater values of the slope wall; recalling the assumption of small slope built into the theory. This situation is observed too if we compare the snap-off time (${\mathsf{\tau }}_s$) of the cases of Figure 8 with those reported of CFD experiments (${\mathsf{\tau }}_{CFD}$) (Figure 3 in [14]); the ratios ${\mathsf{\tau }}_s/{\mathsf{\tau }}_{CFD}$ for ${\mathsf{\alpha }}_{a,b,c}=0.45,0.225,0.0695$ are (a) $2.70$ (b) $2.08$ (c) $1.13$, respectively.

4. Conclusions

The snap-off mechanism in circular capillaries was studied here from a theoretical approach. Particularly, we focused on the models reported that are based on the small-slope approach [9,10] for gas–liquid elastic, gas–liquid and liquid–liquid interfaces, mainly: GR, HE and BD models [9,10,16,17]. Although the rheology of each interface implies important considerations on the formulation of the models, their structures are similar: evolution equations described by highly nonlinear partial differential equations of fourth order in space and first order in time (Equations (8)–(12)). According to the reviewed literature, different mathematical techniques have been employed to solve such models, from the Forward Euler simple method to the use of powerful solvers implemented in commercial software.

We thus present a unifying framework to solve the group of referred models. This framework is based on Fourier pseudo-spectral (FPS) differentiation method and uses the Fast Fourier Transform (FFT) and the inverse FFT (IFFT), as is schematized in Figure 3. To demonstrate the effectiveness of our framework we evaluate some of the most representative model cases studied here and compared the results that we obtained with those reported in the literature. As it can be seen in Figure 4, Figure 5, Figure 6, Figure 7, Figure 8 and Figure 9 there is an acceptable agreement between the results of the snap-off dynamics (shape of the interface profile and snap-off time). The fact that the FPS differentiation method turns out to be suitable for the solution of the models studied here, may be related with the periodic domain of the problem, which results adequate to implement the trigonometric interpolant used in the pseudo-spectral method.

Finally, it is important to mention that despite the fact that the snap-off phenomenon was studied here by mathematical modelling, most of the exposed models (except the HAC model) have been validated experimentally [9,10,13,27,28]. Furthermore, it should be noted that a series of simplifications are made on the models’ formulations, so future research might focus on overcoming these limitations.