Multiphase CFD modeling of front propagation in a Hele-Shaw cell featuring a localized constriction

We study a liquid-gas front propagation in a modulated Hele-Shaw cell by means of multiphase computational ﬂuid mechanics based on the three-dimensional Navier-Stokes equations. In the simulations an obstacle that partially ﬁlls the gap is placed at the center of the cell, and the liquid-gas interface is driven at a constant velocity. We study the morphological differences between imbibition and drainage for a wide range of capillary numbers, and explore how the wetting properties of the constriction affect the amount of liquid that remains trapped in the draining process. We observe increasing remaining volumes with increasing capillary number and decreasing contact angle. The present CFD implementation for a single mesa defect provides insight into a wide number of practical applications.


I. INTRODUCTION
Two-phase fluid displacements in porous media have been extensively studied [1,2] to gain insight into a broad variety of applications that involve interfacial problems, e.g., extraction, filtering, chemical reaction, mass-transfer, and separation. The simplest system involves two immiscible fluids, characterized by their viscosity and surface tension at their mutual interface. A remarkable feature in such systems is the hysteresis appearing during repeated series of drainage-imbibition processes. This is usually evidenced by pressure-saturation trajectories, which follow different paths depending on the flow direction. Since this behavior does not relate to disorder [3], studying a single mesa defect provides new insight for the collective, macroscale hysteresis [4].
In the confined geometry of a Hele-Shaw cell (two closely spaced parallel plates) the interfacial dynamics depends both on the relative ability of the fluids to wet the walls of the cell and on their viscosity contrast. We focus on the simplest case, where a viscous wetting fluid displaces a second less-wetting and less-viscous fluid. In this imbibition process the viscosity contrast stabilizes the interface and makes it effectively stiffer at increasing front velocities; in the opposite displacement direction (drainage) the viscosity contrast makes the interface increasingly softer with velocity and eventually unstable against the Saffman-Taylor viscous fingering instability. This simplest case is essentially governed by the competition of viscous and capillary forces, as represented by a dimensionless capillary number, Ca.
Roughness and surface heterogeneities have an important effect on liquid spreading. In the imperfect Hele-Shaw cell introduced by de Gennes [4], defects of either geometrical or chemical origin deform the fluid-fluid interface and impact the interfacial dynamics. Several studies suggest that the nonlocal nature of the interfacial dynamics and the presence of defects at the interface can be employed to control interfacial stability [5][6][7][8][9]. While keeping the same fluid pair, the instability is delayed by the effect of wetting and increasing the wettability of the invading fluid stabilizes the invasion [9][10][11]. Similarly, gap-thickness modulations and localized gap-thickness defects (constrictions and expansions) have been shown to produce controlled shape distortions of stable fronts [12], as well as front pinning and front depinning (accompanied by abrupt capillary jumps) across sharp-edged defects [3].
Beyond the quasistatic limit corresponding to Ca = 0, the interfacial dynamics in Hele-Shaw cells with gap thickness fluctuations is nonlocal: it is governed by a stochastic integro-differential equation for the gap-averaged front height, with different noise terms accounting respectively for localized changes in capillary pressure, cell permeability, and local mass conservation, all of them arising from the localized changes in gap thickness [13,14]. In this framework, addressing large interfacial deformations and interface pinch-off seem to fall beyond the possibilities of most analytical approaches. Instead, numerical approaches to the complete three-dimensional problem look more suitable.
In this work we study liquid-gas front propagation in a modulated Hele-Shaw cell by means of multiphase computational fluid mechanics based on the full three-dimensional multiphase Navier-Stokes equations. We consider the scenario of a single mesa-shaped defect which modifies locally the pressure and wettability. Our goal is characterizing the morphological differences between imbibition and drainage for a wide range of capillary numbers, and exploring how the wetting properties of a localized constriction influence the amount of liquid trapped in the course of drainage. Our implementations can easily be extended to disordered media, where single defects interact and collectively generate the pressure-saturation relationship.
The simulations provided here demonstrate the capabilities of the current multiphase fluid dynamics methods for simulating the interfacial dynamics of fluid front propagation in quasi-twodimensional systems. Furthermore, we are able to examine the influence of the physical properties of fluid arbitrarily, and alter the wettability of the cell, which has only recently been studied experimentally and has not been fully explored [10,11]. The simulations provide us also with the velocity field, which reveals the generation of streamline vortices during pinch-off.
The remainder of the paper is organized as follows. The theoretical framework and the numerical procedure are presented in Sec. II. The numerical model is compared with well-known experimental results in Sec. III. We perform a systematic study of the impact of wetting properties of fluid-fluid invasion and the capillary number in Sec. IV. Discussions are presented in Sec. V and VI, and the main conclusions drawn in Sec. VII.

II. THEORETICAL FRAMEWORK AND NUMERICAL PROCEDURE
We use a square Hele-Shaw cell of size L x × L y = 63 × 63 mm and gap thickness L z = 0.3 mm, as shown in Fig. 1. The gap modulation consists of a square patch of size l x × l y = 3 × 3 mm and height l z = 0.06 mm placed at the middle of the bottom plate, so that the gap thickness on top of the modulation takes the value h = L z − l z . We consider a two-phase displacement, in which the initial interface is flat and the front velocity imposed at the inlet is V ∞ = (0, V ∞ , 0). First, a viscous fluid with dynamic viscosity μ invades the Hele-Shaw cell at constant velocity through the inlet located at one side, simultaneously displacing the air inside the cell [ Fig. 1(a)]. Once the front has propagated sufficiently far past the obstacle, we reverse the direction of the flow, draining the viscous fluid from the cell at a constant velocity V ∞ = (0, −V ∞ , 0) [ Fig. 1(b)]. The cell is tilted at an angle 0 • β 90 • with respect to the direction of gravity g [ Fig. 1(c)].
We perform full 3D numerical simulations of the physical system in Fig. 1 solving the multiphase Navier-Stokes equations employing the volume of fluid approach. The problem is implemented 084305-2 in the OpenFOAM software package [15] and solved using interFoam [16], extensively used for multiphase simulations [16][17][18]. This numerical solver is suitable to address large interfacial deformation and topology changes of the fluid-fluid interface. The governing flow equations consist of the mass balance equations and the full Navier-Stokes equations including the interface force term f γ defined as where α is the liquid fraction, γ the surface tension. and κ the local interfacial curvature determined by Wettability effects are included in the model through the contact angle θ , where n w and t w are the unit vector normal and tangent to the wall, respectively [17]. The density ρ and viscosity μ fields are defined as where (ρ l , μ l ) and (ρ g , μ g ) are constants in each phase (liquid and gas). The two fluids are considered to be Newtonian, immiscible, incompressible, and isothermal. In the volume-of-fluid method, the liquid fraction α varies smoothly between 0 and 1, and from mass conservation and Eq. (4), it evolves according to ∂α ∂t 084305-3 In interFoam solver, the phase is limited using the Multidimensional Universal Limiter with Explicit Solution (MULES) and, due to the presence of sharp gradients of α in the interfacial region, an additional convective term is added in Eq. (6) to mitigate the effects of numerical diffusion [16]. The momentum equation is solved via a Pressure Implicit with Splitting of Operators (PISO) iteration procedure with nCorrectors (the number of times the algorithm solves the pressure and momentum corrector at each step) equal to 4. The time derivative is solved by an Euler scheme, momentum and mass flux by the Gauss upwind scheme. The system of differential equations is solved using the PCG linear solver with DIC preconditioning, with an absolute tolerance of 10 −10 . In order to keep sharp interfaces, the cAlpha parameter controlling the interface compression is set to 1. Finally, the parameter nAlphaCorr specifying the number of times the alpha field is solved within a time step here equals 3.
The computational domain consists of half Hele-Shaw cell (see Fig. 1) with symmetric boundary condition in the x-axis. The rest of the boundary conditions are no-slip at all walls, constant velocity V ∞ at the inlet, and zero pressure at the outlet. A constant contact angle is used as wetting boundary condition for either θ t or θ b , and the wettabilities of the cell and obstacle are defined in terms of the contact angles θ h and θ o , respectively. Values θ = 0 • (θ = 90 • ) correspond to perfect wetting (superhydrophobic behavior) of the invading fluid, respectively. The rest of the simulated liquid parameters vary according to Table I, which are the most common in practical applications. Under these conditions the capillary number is 2.8 × 10 −6 Ca 4.8 × 10 −4 and the Reynolds number is 0.0012 Re 0.14.
After a convergence study, we find that using 158 760 mesh elements allows us to reach accurate solutions and reasonable computation times. It is worth pointing out that all the implicit interface capturing methods are susceptible to generation of unphysical flow at the interface only due to numerical reasons [21]. However, Harvie et al. [22] have shown that the generation of spurious currents is only secondary for moving interface, even though the following time step criterion is applied [16,23]: where are the viscous and capillary time scales, respectively, and x the cell size.

III. QUALITATIVE COMPARISON WITH EXPERIMENTS
In this section we compare the efficiency of the solver with theoretical and experimental results on the well-known Saffman-Taylor instability. We consider the full-wetting case in which both contact angles θ h = θ o = 0 • and the fluid-fluid pair is oil 1-air. Figure 2(a) shows the position of the interface at different times during imbibition (left) and drainage (right) displacements for the velocity V ∞ = 0.2 mm/s (capillary number Ca = μV ∞ /γ 5 × 10 −4 ) and β = 0 • . During the imbibition stage the cell gets filled uniformly by oil, the interface remains flat and it is only slightly modulated during the transition over the obstacle. Furthermore, the final interface configuration shows no signs of the presence of the obstacle, nicely following recent experimental observations [3]. On the contrary, two distinct aspects manifest during drainage. A Saffman-Taylor viscous fingering instability takes place, as known to occur when a more viscous fluid is displaced by a less viscous one in a horizontal cell. In addition, the interface pins at a point P 1 (see Fig. 1) on the obstacle and at the cell walls, whereas the rest of the interface continues retreating. As a result, some of the liquid remains on top of the obstacle, and the interface is squeezed by the surrounding air until it pinches off.  . This leaves the interface pinning at P 1 unchanged, while the fingering instability is more pronounced at higher driving velocity, as expected. In this case the interface is not squeezed around the obstacle and breaks up only when parts of it reach the inlet.
The amount of fluid remaining on top of the obstacle V changes with the inlet velocity V ∞ . A thin coating layer of liquid remains on the cell plates, whose thickness increases with V ∞ . The presence of thin film layers covering the cell walls during the drainage stage has been experimentally observed by Liu and coauthors [24].
Assuming that the obstacle acts only as a local perturbation to the interface, we may confirm that the instability observed here is indeed viscous fingering by comparing the observed number of fingers to that expected to appear during drainage, computed by following Park and Homsy [25,26]. By applying linear stability analysis to the flat interface, the most unstable mode has a dimensionless wave number given byk Now, when V ∞ increases, the distance between fingers decreases. This prediction aligns with the results of our simulations, shown in Fig. 2.
Tilting the Hele-Shaw cell an angle β adds an effective gravity force. This force can be large enough to suppress the viscous fingering instability, and then all the liquid remaining on the sides of the cell is drained. In addition, the amount of liquid remaining on top of the obstacle decreases (see Fig. 3), as one might expect due to the additional body force moving the fluid. We can estimate the critical angle β c for which the fingering disappears as [3,14] β c = arcsin 12 where μ = μ l − μ g and ρ = ρ l − ρ g . For a Hele-Shaw cell of gap thickness h = 0.46 mm, the critical angle is β c 3 • , which is smaller than the angle used by Planet et al. [3] in their experiments to prevent viscous finger formation in drainage.

IV. PARTIALLY WETTING OBSTACLE
We consider now the scenario where the obstacle is only partially wetting, θ o = 0, and the cell remains completely wetting, θ h = 0. Now the pressure around the obstacle slightly changes according to Eq. (2) and, depending on the value of θ o , the interface may or may not be pinned on it. In the case of oil 1-air and θ o > 20 • , the interface experiences smooth corrugations at the obstacle, but the front does not get pinned. Figure 4 shows a comparison of two cases, for θ o = 10 • and θ o = 30 • . Note that the perturbation at the interface is stronger than in the fully wetting case.
Replacing the oil with water makes the interface more curved, but the amount of liquid that remains trapped on top of the obstacle during the liquid retraction (drainage) still depends on the contact angle θ o . In Fig. 5 we show snapshots of the simulation at different times during the injection (left column) and withdrawal (right column) of water in the cell for selected cases having different contact angles (a) θ o = 45 • and (b) θ o = 60 • . As the contact angle θ o increases, the volume of liquid remaining on top of the obstacle V decreases until it reaches zero. When that occurs, the shape of the interface in the retraction phase (drainage) behaves similarly to the injection phase (imbibition), where after surpassing the obstacle the interface is no longer affected by it. Furthermore, we note that the injection phase is less influenced by the change in the contact angle, complying with the fact that the interface is stiffer in that particular direction of flow.  Figure 6 shows the corresponding evolution of the front for (a) ether and (b) ethanol with θ o = 45 • . As soon as the interface makes contact with the obstacle during drainage, at point P 1 (see Fig. 1), the pinning force on the interface produced by the obstacle competes with the weakly restoring force of the interface, giving rise to the break-up of the interface below the obstacle close to point P 2 .
Interestingly, when the capillary number increases significantly the behavior of the interface changes completely. Figure 7 shows the interface evolution for an oil-air interface with the oil contact angles set to (a) θ = 45 • and (b) θ = 60 • . During imbibition, the interface behavior is independent of the contact angle θ o . The formation of viscous fingers in drainage, however, makes the velocity of the fluid at the lateral walls to slow down abruptly and the interface to pin at point P 1 of the obstacle. The rest of the interface keeps advancing and begins to curve, but in contrast with the previous liquids there is no pinch-off until the interface finally reaches the inlet boundary. At that moment the break-up takes place and the liquid starts to relax slowly over the obstacle until it reaches a symmetric shape just as in the cases discussed earlier [ Fig. 7(a)]. When the contact angle θ o is larger, the interface does not pin at the obstacle, but instead carries on moving at a velocity different than the rest of the interface. Finally, when the interface reaches the boundary it breaks up and a droplet of liquid remains inside the cell, but it is not placed over the obstacle, as shown in Fig. 7(b).

V. PINCH-OFF AND REMAINING VOLUME
The simulations show that the pinch-off process depends on the wettability conditions, gap modulation and capillary number. In addition, we see that before the breaking of the interface the fluid is squeezed out from the neck, while the interface retracts. The condition for pinch-off to take place is that the time invested by the interface at the obstacle to travel before the neck collapses is larger than the time required for the neck to pinch-off. If we consider a thin film of thickness x n (neck width), it is easy to see that the pinch-off time is τ p ∼ μx n /γ .
The amount of volume that remains on top of the obstacle V is related with the neck width x n and slightly changes with the contact angle θ o . To include the different liquids used here, we define the capillary number and measure the volume V and neck width x n . Figure 8 summarizes how the dimensionless volume V/(l x x 2 n ) changes with the obstacle wettability properties through the capillary number defined by Eq. (12).

VI. TRANSITION FROM NO BREAKUP TO BREAKUP
The imbibition of a high-viscosity fluid in a Hele-Shaw cell filled with a less viscous fluid shows that the interface is morphologically stable, and no significant changes modify its shape. More interesting is the drainage displacement, where eventually the retracting liquid front remains pinned on the obstacle until pinch-off occurs. Further insight about this process can be gained from Fig. 9(a), which shows a simulation image of the necking process leading to interface break-up in the case of retraction of ethanol. The displayed velocity field shows the fluid squeezing out from the neck to the obstacle, while the interface retracts. Once the interface reaches the obstacle, a recirculation zone is created. The interface pushes the recirculation zone along the cell as it advances towards the inlet. In contrast to the case of ethanol, Fig. 9 zone for retracting oil increases. Due to this, the width of the neck close to the obstacle does not shrink, and therefore no pinch-off is observed.
To quantify the effect, we use the vorticity defined by where is the volume occupied by the liquid. Figure 10 shows the vorticity for several cases as a function of time. At early times, the vorticities behave smoothly, while at later times they present huge variations. The transition between these two regimes corresponds to the instance when the interface reaches the obstacle, i.e., the moment when the recirculation zone is generated (Fig. 9). Two different aspects are interesting in the previous figure. First, the absolute value of the vorticity increases with capillary number, which agrees with the pinch-off volume as described above, whereas the contact angle at the obstacle does not modify the vorticity. Second, a significant change appears when we move from a system where pinch-off does not occur to one where pinch-off takes place around the obstacle. For instance, the vorticity for oil-air and θ h = 0 (Fig. 7) is one order of magnitude bigger than the case with θ h = 0 (Fig. 2).

VII. CONCLUSIONS
The present study considered the imbibition and drainage of liquid in a Hele-Shaw featuring a single constriction, analyzed using multiphase computational fluid mechanics based on the threedimensional Navier-Stokes equations as implemented in the interFoam solver of the OpenFOAM library. We have explored the dependence of wettability conditions and the capillary number are a useful tool to reproduce nontrivial behaviors arising from the interplay between viscous and capillary forces.
During the first stage, the cell is filled by the more viscous fluid, which displaces air. We have shown that during this imbibition stage the interface is stiff and no deformations are observed. However, exception to this is observed at large obstacle contact angles, where complete imbibition is no longer reached due to the obstacle-liquid repulsion.
An individual obstacle seems to not greatly affect the front, i.e., the front keeps straight after the interface overtakes it. The permeability effect due to the gap thickness variation is well compensated by the locally increased pressure difference atop of the obstacle, which pushes the interface ahead. Therefore the reverse process, the drainage stage, reveals the true effect of the obstacle in the hysteresis behavior.
During the drainage stage, the model reproduces the well-known undulations of the interface. The fingering instability and the presence of the obstacle are decoupled: while the amplitude of the instability increases as V ∞ increases, the interface is always pinned at the obstacle. Additionally, experiments have shown the suppression of the fingering instability due to gravity, and that these undulations completely disappear at tilting angles above a critical one. We found that our simulations including the effects of gravity reproduce this behavior very closely. Even the critical tilting angle observed was rather close to its experimental predicted value.
The simulations allowed for a detailed analysis of the flow field. We observed that the presence of the obstacle creates a recirculation of the velocity field. While the liquid is retracting, the amount of recirculation around the obstacle is reduced and the recirculation zone is pushed to a small region. Then, two different situations may arise: either the interface breaks up close to the obstacle, or it breaks up at the inlet. Both situations depend on the capillary number and wettability conditions.
To conclude, an analysis of the stability of the front based on the restoring forces provided by surface tension and gravity, the destabilizing force of the viscous pressure gradient, and the geometry of the obstacle reveals the wettability conditions that give rise to pinch-off in the drainage displacement. Therefore, the results shown here can be used to either avoid or control the motion and stability of liquid fronts in microscale geometries, and may be useful also to understand more complicated situations. Moreover the stability conditions can be used to design new platforms in microfluidics for the manipulation of tiny liquid volumes, with tailored spatial covering.

084305-12
Imbibition of yield stress and other complex fluids is an emerging topic in soft matter partially due to its practical applications in oil recovery. This simulation method allows for a natural extension of the tools to simulate non-Newtonian fluid front propagation in a Hele-Shaw cell and in porous random media.