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

: The snap-o ﬀ is an instability phenomenon that takes place during the immiscible two-phase ﬂow in porous media due to competing forces acting on the ﬂuid phases and at the interface between them. Di ﬀ erent theoretical approaches have been proposed for the development of mathematical models that describe the dynamics of a ﬂuid / ﬂuid interface in order to analyze the snap-o ﬀ 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 ﬂow at pore scale in circular capillaries for pure and complex interfaces. The models consist of evolution equations; highly nonlinear partial di ﬀ erential equations of fourth order in space and ﬁrst order in time. Although the structure of the models for each type of interface is similar, di ﬀ erent 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 di ﬀ erentiation 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-o ﬀ .


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].

Figure 1.
Snap-off phenomenon illustration in a water wet porous media which contains residual oil ganglia (in black). Snap-off occurs when: (a) the front of the non-wetting phase flows near the pore throat; (b) Then the wetting phase is displaced by the non-wetting one, leaving a thin film at the wall; (c) A liquid collar of the wetting phase grows due to the capillary pressure difference at the crest and the capillary throat; (d) The collar breaks up the displacing phase due to the surface tension and a separate ganglion is formed. 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 Figure 1. Snap-off phenomenon illustration in a water wet porous media which contains residual oil ganglia (in black). Snap-off occurs when: (a) the front of the non-wetting phase flows near the pore throat; (b) Then the wetting phase is displaced by the non-wetting one, leaving a thin film at the wall; (c) A liquid collar of the wetting phase grows due to the capillary pressure difference at the crest and the capillary throat; (d) The collar breaks up the displacing phase due to the surface tension and a separate ganglion is formed.

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 nonwetting 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 twophase flow at pore-scale. Although any smooth function can be used to represent this particular geometry, a sinusoidal function is typically used [17]: where  is the radial coordinate of the capillary wall, x refers to the axial coordinate, a is the ratio between the throat ( t R ) and capillary ( c R ) radii and  is a geometric parameter which represents the ratio between the capillary radius ( c R ) and half of the capillary length (), ( 2) c R  .

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]: Figure 2. Geometric diagram that shows circular constricted capillaries and the configuration of the liquid phases. The main geometric elements are the radius of the capillary R c ; the throat radius R t , and the length of the capillary . The configuration of the phases consists of a core-annular flow, in which the external phase (with subscript 2), with constant thickness δ and viscosity µ 2 , wet the wall of a pore described by λ(x); and the internal phase, with viscosity µ 1 , is surrounded by the external phase, with which it shares the interface κ(x, t). The volumetric flow Q, moves from left to right. 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]: where λ 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 α is a geometric parameter which represents the ratio between the capillary radius (R c ) and half of the capillary length ( ), α = R c /( /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]: where P is the pressure, subscripts 1 and 2 denote, respectively, the non-wetting and wetting phase, µ 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, κ the radius of interface and σ the surface tension.
The drop breakup problem needs to determine the interface position change over time ∂κ(x, t)/∂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: 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 κ. 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.

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.

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 κ, 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 = λ → u 2 = 0 ) and balance of forces at the interface (r = κ → µ 1 du 1 /dr = µ 2 du 2 /dr) 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 (µ 1 = 0), shear stress of the core is neglected µ 1 du 1 /dr = 0; consequently r = κ → µ 2 du 2 /dr = 0 . In Hoyer and coworkers formulation [16] for elastic interfaces, surfaces stress (T = σ + Kξ) in axial (xx) and transversal (θθ) directions are included as function of the surface tension σ, the dilational interfacial elasticity K and the corresponding surface strains ξ xx and ξ θθ (defined by geometric relationships). Note that for inelastic cases (K = 0), the surface stress becomes T xx = T θθ = σ. At the interface the liquid film shear stress is balanced by the net surfaces stress in axial direction, then: r = κ → µ 2 du 1 /dr = dT 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 = σ 1/κ + ∂ 2 κ/∂x 2 ; while for an elastic interface the expression is For fluid arbitrary viscosities (µ 1 0), the velocity profile of the core is equal to that of the annular film at the interface ( r = κ → u 1 = u 2 ); additionally no variation of the radial velocity derivative on the axes is postulated (r = 0 → du 1 /dr = 0). 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: The nondimensionalizing is performed according with: where * denotes dimensionless variables and δ is the film thickness. For the viscous core case, the dimensionless equations use the viscosity of the non-wetting phase µ 1 instead of the wetting phase µ 2 , and the dimensionless time is τ = t/(µ 1 R c /σ).

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]: if K * = 0, i.e., for inelastic interface, it then results in an extension of the GR model [16]: The Gauglitz and Radke model (GR) for the gas-liquid interface is [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., λ * (x * ) = 1, the GR model results in: The Beresnev and Deng model (BD) for liquid-liquid interfaces is [17]:

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 δ 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.

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 (Ca) : 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:

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: 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: Mathematics 2020, 8, 1941 8 of 15 Periodic boundary conditions for the BD model are imposed, that is, the interface radius is equal at the boundaries for all simulation time: 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.

Numerical Procedures
The models presented in Sections 2.3.2 and 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 * = [−L/2 L/2] divided on N nodes spaced for a distance h = 2π/N, period T = 2π, wave numbers k entire multiples of π/h = N/2, is stablished by: for all k = −N/2 + 1, −N/2, . . . , N/2; whereκ 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: 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 τ, κ * = ∂κ(x * , τ)/∂τ, on the domain x * = [−L/2 L/2], begins with the compute of the DFT by the FFT (κ * ). Subsequently the n-order derivative ofκ * is performed on the Fourier space by using the associated wave numbers k. Finally, the derivative κ * is obtained via the IFFT.
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  Once obtained the n-order spatial derivatives of the interface position *  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.

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 0.2 a  , 0.1  . The initial condition is defined by Figure 4 shows the evolution of the interface radius *  Once obtained the n-order spatial derivatives of the interface position κ * 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.

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, α = 0.1. The initial condition is defined by Ca = 0.1. First, Figure 4 shows the evolution of the interface radius κ * 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 × 10 −4 . In the case of a straight capillary, the GR model is evaluated for a dimensionless initial film thickness ε = 0.12 and a dimensionless wavelength of fastest growing disturbance of λ * max = 2 3/2 π. 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. ). 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 ). While in Figure 5, it is related to the interface profiles moments before the end of the simulation for different values of * K .  ). 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 methods for the HAC model (gas-liquid elastic interface). The results correspond to case: ). 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 ). 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 ) and FPS ( ). 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 ) methods where a = 0.2, α = 0.1, and K * = 0 ). 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 ). 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   ). 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 same geometry used on the previous case but with initial conditions defined by . In the case of a straight capillary, the GR model is evaluated for a dimensionless initial film thickness 0.12  and a dimensionless wavelength of fastest growing disturbance of 32 * max 2    . 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.   In the case of a straight capillary, the GR model is evaluated for a dimensionless initial film thickness 0.12  and a dimensionless wavelength of fastest growing disturbance of 32 * max 2    . 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.   .
) and FPS ( the same geometry used on the previous case but with initial conditions defined by 8.94 10 Ca x  . In the case of a straight capillary, the GR model is evaluated for a dimensionless initial film thickness 0.12  and a dimensionless wavelength of fastest growing disturbance of 32 * max 2    . 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.  In the case of a straight capillary, the GR model is evaluated for a dimensionless initial film thickness 0.12  and a dimensionless wavelength of fastest growing disturbance of 32 * max 2    . 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.  ; τ = 3000 the same geometry used on the previous case but with initial conditions defined by 8.94 10 Ca x  . In the case of a straight capillary, the GR model is evaluated for a dimensionless initial film thickness 0.12  and a dimensionless wavelength of fastest growing disturbance of 32 * max 2    . 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.  ; τ = 5300 In the case of a straight capillary, the GR m 0.12  and a dimensionless wavelength o we compared the evolution of the interfac Galerkin FE method, to those obtained with the results of this work applying the FPS m  In the case of a straight capillary, the GR model is evaluated for a dimensionless initial film thickness 0.12  and a dimensionless wavelength of fastest growing disturbance of 32 * max 2    . 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.   In the case of a straight capillary, the GR model is evaluated for a dimensionless initial film thickness 0.12  and a dimensionless wavelength of fastest growing disturbance of 32 * max 2    . 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.  ) and its comparison with those reported by Gauglitz and Radke (1988) with Galerkin-FE ( In the case of a straight capillary, the GR model is evaluated for a dimensionless initial film thickness 0.12  and a dimensionless wavelength of fastest growing disturbance of 32 * max 2    . 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.  In the case of a straight capillary, the GR model is evaluated for a dimensionless initial film thickness 0.12  and a dimensionless wavelength of fastest growing disturbance of 32 * max 2    . 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.  ), for ε = 0.12 and λ max = 2 3/2 π. The profiles are shown for: τ = 885.4 In the case of a straight capillary, the GR model is evaluated for a dimensionless initial film thickness 0.12  and a dimensionless wavelength of fastest growing disturbance of 32 * max 2    . 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.  ; τ = 2766.2 In the case of a straight capillary, the GR model is evaluated for a dimensionless initial film thickness 0.12  and a dimensionless wavelength of fastest growing disturbance of 32 * max 2    . 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.  ; τ = 12870.4 In the case of a straight capillary, the GR model is evaluated for a dimensionless initial film thickness 0.12  and a dimensionless wavelength of fastest growing disturbance of 32 * max 2    . 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, δ * = 0.5R t * , µ 1 /µ 2 = 10 and different values of α. 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 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   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 ( 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: 12   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   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 ( 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: 12   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   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 ( 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: 12 0.1    . The rest of the parameters are those of the base case, but with 5 3 10 Ca x   and 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 s ) methods. The results correspond to the simulations performed with Ca = 10 −3 , a = 0.6, δ * = 0.5R t * , µ 1 /µ 2 = 10 and the geometric slopes: α = 0.45 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   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 ( 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: 12 0.1    . The rest of the parameters are those of the base case, but with 5 3 10 Ca x   and 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 s ; α = 0.225 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   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 ( 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: 12   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 (α = 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: µ 1 /µ 2 = 0.1. The rest of the parameters are those of the base case, but with Ca = 3 × 10 −5 and α = 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 τ s ; since this is increased from τ s = 39.54 to τ s = 600.16. Beresnev and Deng (Figure 3c and Figure 7 in [14]) report for these cases τ sBD = 38.0, 585.0, respectively.
As it may be observed on the series of Figures 4-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 (Figures 5 and 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 (Figures 8 and 9). This can be explained by the high non-linearity of the BD model (Equation (12)) and the periodic boundary conditions imposed (Equation (19)).  (Figures 5 and 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 (Figures 8 and 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 1  , so the smaller  , the better the approximation. On Beresnev et al. (Table III in

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. As it may be observed on the series of Figures 4-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 (Figures 5 and 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 (Figures 8 and 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 1  , so the smaller  , the better the approximation. On Beresnev et al. (Table III in

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. As it may be observed on the series of Figures 4-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 (Figures 5 and 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 (Figures 8 and 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 1  , so the smaller  , the better the approximation. On Beresnev et al. (Table III in

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. As it may be observed on the series of Figures 4-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 (Figures 5 and 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 (Figures 8 and 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 1  , so the smaller  , the better the approximation. On Beresnev et al. (Table III in

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. As it may be observed on the series of Figures 4-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 (Figures 5 and 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 (Figures 8 and 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 1  , so the smaller  , the better the approximation. On Beresnev et al. (Table III in

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.
It is important to remark that the small-slope based models assume α 1, so the smaller α, 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 (τ s ) of the cases of Figure 8 with those reported of CFD experiments (τ CFD ) (Figure 3 in [14]); the ratios τ s /τ CFD for α a,b,c = 0.45, 0.225, 0.0695 are (a) 2.70 (b) 2.08 (c) 1.13, respectively.

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 Figures 4-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.
Funding: This research was funded by the National Council for Science and Technology (CONACYT) in Mexico, grant number 25345.

Conflicts of Interest:
The authors declare no conflict of interest.