CONTROL OF FINGERING INSTABILITY BY VIBRATIONS

In applications involving the injection of a fluid in a porous medium to displace another fluid, a main objective is the maximization of the displacement efficiency. Displacement fronts moving in porous media are subjected to hydrodynamic instability when a liquid of low viscosity displaces a high-viscosity liquid and consequently finger-like structure forms along the interface. This finger instability is usually undesirable in technical applications and natural filtration processes. We discuss the external periodic forcing as one of the promising ways to control the instability and perform numerical simulation of an initially spherical drop in a porous media under vertical vibrations. The drop is a favorable object to study since in this case one can observe the effect of vibrations on fluid interface domains inclined by different angles with respect to vibration axis. It is shown that under vibrations small-scale perturbations of interface are suppressed and in the case of vibrations of large enough intensity the drop becomes stable. The stability criterion is derived. Mathematics Subject Classification. 76E17, 76S05. Received November 4, 2020. Accepted May 19, 2021.


Introduction
In the framework of the conventional Darcy model a planar displacement front in a porous medium was found to be unstable if the displacing fluid has lower viscosity than the displaced one [2]. The instability is associated with the development of penetrating fingers of the less viscous fluid. They grow due to the lower friction with porous matrix. The process is virtually unimpeded by viscous momentum transfer, which is negligible in porous media. This mechanism explains the short-wave nature of the instability.
This finger pattern is observed in different technical applications and natural filtration processes, for instance oil and gas production [2,9,20], geological CO 2 storage [5,7], development of low-temperature polymerelectrolyte fuel cells [24], water seepage [6], and ink spreading [1]. This instability is usually undesirable and in the present paper we analyze the use of vibrations to control the stability of fluid-fluid interface in a porous medium. For example, the displacement of oil by water injected into an oil reservoir is unstable, which ultimately results in the formation of retained oil pillars surrounded by water. The problem of the extraction of the retained oil is a topical one, since the solution of this problem could noticeably increase the recovery efficiency. Another example is the production of titanium [11], according to the technology the producing magnesium drops of molten magnesium have to be extracted from a porous titanium sponge saturated with liquid magnesium chloride.
In practice the displaced fluid often forms compact inclusions surrounded by the displacing fluid. In this case one should consider displacement fronts on both sides of the inclusion. The investigation of the motion of localized liquid inclusions in a porous medium saturated with another liquid under the action of stationary forces is of interest for the theory of hydrodynamic stability of multiphase systems. In most problems of this type, the equations of filtration of two-component liquid are used. However, in the initial stages of displacement the mixing of the liquids can be ignored and the phase interface can be considered as macroscopic. Moreover, an effective phase interface (saturation jump) may also be formed during the evolution of the mixture [2].
In the framework of the linear stability theory the dynamics of the steady settling of spherical inclusion of one liquid in a porous medium saturated by the other liquid was studied in [11]. It was shown that monotonically increasing perturbations are localized at the forefront of moving droplet, in agreement with the Rayleigh-Taylor mechanism. On the other hand, a two-parameter family of eigenmodes was found, with the spectrum of eigenvalues filling the entire complex plane.
The suppression of the finger instability might be useful to improve the efficiency and control of a number of physical, biological, and technological processes related to the viscous fingering phenomena. The stabilization of fluid-fluid displacements can be achieved using wettability alteration [23] or implementing the porous media of certain structure [18]. Another promising ways to control the instability is the implementation of external periodic forces. Vibration effect on a stability of plane displacement front in a porous medium was studied in [14]. It was found that the vibration stabilizes the plane displacement front with respect to shortwave perturbations that usually have the highest grow rates. Numerical modeling in [8] show that moving drop can be stabilized by the modulation of external pressure gradient.
In the present paper the numerical simulation of liquid drop dynamics in a porous medium under gravity is performed and the effect of modulation of external forcing is analyzed. Calculations are performed using the Level set method. The method was developed in [4,13,17] where interface between fluids of different properties is represented as transition layer of finite thickness such that the fluids are considered as one medium with sharp change of properties in the transition region. In [22] a new enhanced treatment of the Level set method was introduced that allows to perform simulation for fluids having large density and viscosity ratios.
The paper is organized as follows. In Section 2 we formulate the problem in terms of nondimensional governing equations and boundary conditions and introduce the nondimensional parameters. In Section 3 we present the results of calculations of a drop dynamics in gravity field. In Sections 4 and 5 the influence of vibrations is considered analytically and numerically. Section 6 contains the discussion and concluding remarks.

Problem statement
Let us consider an infinite porous medium filled with a fluid of density ρ 1 and viscosity η 1 . In the porous medium, a spherical drop of another fluid with density ρ 2 and viscosity η 2 moves under gravity. The drop on the one hand is commonly encountered object in a porous media, on the other the study of its dynamics allow to investigate the influence of vibrations on the fluid interface inclined by different angle with respect to vibration axis. Certainly the conclusions about displacement front could be made only if perturbations of the drop shape are much smaller than the drop radius. We will consider the influence of vertical vibrations on the drop dynamics because in that case asymmetric treatment is possible. Both fluids are assumed to be incompressible.
We neglect the smearing of the phase interface and consider the processes occurring on time scales smaller than the characteristic time of interpenetration of the fluids through the pores. Thus, following to [11,12] it is assumed that the miscibility of the fluids and the size of the drop are such that the interface can be considered sharp. Furthermore, the drop is sufficiently large in comparison to the typical pore size, which allows us to use the pore-averaged description. It is worth to note that the first and the second fluids can be either pure liquids or homogeneous mixtures. One possible configuration which can be described by the suggested model is the mixture of two liquids saturating the entire porous medium in the case that their saturation experiences a jump across some surface while in all other volume the saturation is homogeneous.
It is well known that in a porous medium, capillarity leads to a diffusive front with a smooth saturation profile as shown in Figure 1, where both phases co-exist on a microscopic scale [3,10]. On the other hand smooth saturation profile will only increase the front stability, as shown for example in [3,21]. Therefore, dealing with front stability analysis we can restrict ourselves to the situation corresponding to the most unstable configuration when the possible effect of capillary forces is neglected.
The dynamics of a drop without vibrations is described by Darcy model and the continuity equation in each liquid [12,14,16]: where p j is the pressure, v j is the filtration velocity defined as the velocity averaged over the volume occupied by fluid, K and are the permeability and porosity coefficients, η j is the dynamic viscosity, ρ j is the liquid density, g is the acceleration of gravity, i is the unit vector directed downwards, the subscript j denotes inner "1" or outer "2" fluid. Under vertical vibrations additional periodic force appears in Darcy equation (2.1), in addition, the inertial term is needed to be taken into account to describe correctly the effect of vibrations of high frequency. Thus, the governing equations describing the dynamics of a drop under vibrations read: here a and ω are the amplitude and the angular frequency of vibrations. At the interface the kinematic condition and the condition of continuity of the normal component of the velocity are imposed [11,12,14]:  The assumptions made in this section are required to obtain the analytical analysis of a front and a drop stability presented in Sections 4 and 5. The other part of the paper is devoted to numerical simulation of drop dynamics. To simplify the comparison of analytical and numerical results the calculations were carried out using the same assumptions as for the theoretical analysis. In particular, the sharp interface approach was used.

Numerical method
At the first stage we performed numerical simulations of the drop dynamics in gravity field in the absence of vibrations. A steady-state regime of gravity-induced drop settlement was studied in [11] where in the framework of linear stability analysis it was shown that the monotonously growing perturbations are localized at the leading front of the moving drop and the oscillatory perturbations have the form of waves traveling along the surface.
Obviously numerical modeling should be performed in the domain of finite size. Therefore the calculations were carried out in a cylindrical cavity of height H and radius R bounded by the vertical rigid wall (see Fig. 2). At upper and bottom boundaries the conditions of zero normal derivative of velocity were set while at rigid boundary the condition of impermeability is imposed. We suppose here that the droplet is small enough, so the influence of boundaries on the drop dynamics (and especially its stability) is small. The choice of rigid wall boundary condition is caused by quick numerical convergence.
Numerical modeling of the liquid droplet motion is performed using the Level set method [22]. The basic idea of the method is the introduction of the distance function -an auxiliary function, which has different signs in different phases, and zero value of this function corresponds to the interface location. At the same time the interface is treated as the transition layer with sharply varying parameters. As a result, a two-phase system is described as a single medium with the parameters (density, viscosity) depending on the distance function.
In accordance with the idea of the Level Set method, the governing equations should be supplemented with the transport equation for the distance function The fluid density and dynamic viscosity are calculated using where H (φ) is Heaviside function defined by the expression where is half-thickness of the transition layer. We assumed that = n∆x, ∆x is the distance between the nodes of computational mesh, n = 3 is the number of nodes located in the half-thickness of the transition layer.
As the scales of the problem, we chose R for the length, Kg (ρ 2 − ρ 1 ) / (η 2 + 0.5η 1 ) for the velocity, ρ 2 gR for the pressure, ρ 2 for the density and η 2 for the viscosity. Note, that the chosen scale of the velocity corresponds to known analytical formula for the velocity of steady motion of spherical droplet in porous medium under vibrations [11]. Performing the standard nondimensionalization procedure of the equations (2.1)-(2.2), (2.5)-(3.4) we obtain The problem contains the following dimensionless parameters: λ = ρ 1 /ρ 2 is the ratio of liquid densities, µ = η 1 /η 2 is the ratio of viscosities, A = µ + 0.5/(λ − 1) is the dimensionless viscosity, ς = H/R is the height of the cavity, and r d = r/R is the drop radius. The Level Set method requires high accuracy treatment near the interface of the drop where the parameters of medium sharply change. Taking this into account we implemented the adaptive mesh refinement algorithm that allows maintaining fine enough mesh near the interface during all time of calculations. The algorithm is realized using Paramesh libraries [15] that supports parallel computing based on Message-Passing Interface and dynamic load balancing procedure for distributing the workloads across multiple computing resources. Thus, the mesh is automatically changed to minimize diffusion of the grid near the interface (see Fig. 4). In addition, at each time step, the corrections of the position of the interface is performed such that the thickness of the transition layer is kept constant [22]. As the result, the variation of the drop mass during calculations did not exceed 1%.
We introduced the stream function to describe incompressible flow with axial symmetry. Poisson equation for the stream function was solved implicitly by Generalized minimal residual method [19]. The time derivative was approximated by the second-order scheme. An approximation of the spatial derivatives was carried out by finite volume method based on the integral form of the conservation equations.

Calculation results
Calculations were performed for cavity with ς = 4 and the droplet radius r d = 0.4. The initial coordinate of the drop mass center is z d = 3ς/4. In Figure 3  Calculations performed for different viscosity ratios of settling and rising drop show that the droplet is unstable and instability develops at the forefront of moving drops regardless of the ratio of fluids viscosities and densities. The results of numerical simulations are in good agreement with the conclusions made in [11]. However, the amplitude of the oscillatory perturbations forecasted in [11] is apparently small because waves traveling along the surface were not observed in numerical modeling.
It would seem worth to note that the drop is unstable with respect to perturbations with the arbitrarily large wavenumbers and the growth rate of small-scale perturbations is the greatest [2,8,14] therefore the result of numerical simulation depends on the parameters of computational mesh. With the increase of the mesh resolution the perturbations of smaller scale become capable to be resolved. The rate of their growth is larger then growth rate of other perturbations so the evolution of the drop shape changes significantly. This effect is illustrated in Figure 4.

Front stability under vibrations
Effect of linearly polarized vibration on the stability of plane displacement front in a porous medium was studied in [14]. Following this paper let us consider the problem of the stability of the motion of a plane displacement front traveling at a uniform velocity U under the action of harmonic vibrations with amplitude a and frequency ω. The front motion could be caused by pressure gradient, gravity, other volumetric force or just inertia. In the framework of considered problem the reason of the front motion is not important.
The flow of each fluid through the porous medium can be described by Darcy equations with the additional term describing vibrations. Let us combine the origin with the displacement front, the z axis being directed perpendicular to the front and the x axis along it. In a traveling coordinate system moving together with the displacement front the governing equations read here U is the velocity of the front motion. From the equation (4.1) one can see that the vorticity (curl v j ) can only damp with time therefore we can introduce the velocity potential ( v j = ∇φ j ) and the governing equations (4.1), (4.2) take the form Boundary conditions (2.5) should be rewrote as The Laplace equation has the following solutions damped at infinity: Table 1. Physical parameters of a water-oil system saturating a clay soil [2,14]. Here, the dot is the derivative with respect to time. Substituting these equations into governing equations, we obtain the Mathieu equation for B In the dimensionless form it readsB In high-frequency limit, applying the standard averaging procedure to equation (4.9), we obtain the following expression for perturbation growth rate [14] In the dimensional form this expression becomes In the absence of vibrations the second term in expression (4.12) vanishes and the stability of the displacement front depends only on viscosity ratio. One can see that under vibrations short-wave perturbations of the displacement front are suspended. Consequently, when a lower viscosity fluid displaces a fluid of higher viscosity, only shortwave perturbations can be stabilized by vibrations and the critical wavenumber of this vibrational stabilization is In the case of vibrations with finite frequencies the equation (4.9) should be solved numerically. The numerical calculation were carried out for the parameters of a water-oil system saturating a clay soil (see Tab. 1). In Figure 5 the dependence of perturbations growth rate on dimensionless wavenumber is described. The results obtained numerically for finite frequencies are in a good agreement with the formula (4.11) derived using high frequency assumption. Numerical solution shows that vibrations of finite frequency also stabilize the short wave perturbation of the displacement front.
For a drop of small enough size one could expect that the vibrations are able to stabilize the drop completely.

Numerical modeling of a drop dynamics under vibrations
The numerical simulations of drop settling under vertical vibrations was performed using the approach and the algorithm described in the Section 3 for the parameter of oil-water system (see Tab. 1). The influence of vibrations was taken into account by introduction of inertial force in the Darcy equation. The time step was taken more then 100 times smaller than the period of vibrations. Calculations were performed using a fine mesh with numbers of mesh nodes per unit of length near the interface equal 256.
Calculations confirmed the results of linear stability analysis of front displacement. The small-scale perturbations of interface are suppressed by vibrations and only large-scale perturbations develop (see Fig. 6a). With the increase of amplitude or frequency of vibrations the intensification of the stabilizing effect is observed (Fig. 6b). For large enough intensity of vibrations the drop becomes stable (Fig. 6c). It should be noted here that we neglect capillarity in our analysis while for oil-water mixture it is important. However, as is it was mentioned in Section 2 in a porous medium, capillarity increases effective diffusivity of the front as shown in Figure 1 [3, 10] and smooth saturation profile will only increase the front stability [3,21]. Therefore, the stability criteria without capillarity is more strict when it should be for realistic oil-water mixture.
The criterion of the suppression of viscous fingering and the drop stabilization can be obtained using the expression for critical wavenumber (4.13). If the size of perturbations which are suppressed by vibrations is greater than the drop diameter then we can expect stabilization of the drop. Thus, we can state that the drop  diameter should be less than the perturbation wavelength, therefore k cr = 2π 2r = π r (5.1) As it was mention above, the velocity of spherical drop of one fluid surrounded by the other fluid in porous media can be calculated using [11,12] U = Kg (ρ 2 − ρ 1 ) (η 2 + 0.5η 1 ) (5.2) Substituting equations (5.1), (5.2) in (4.13) we obtain the following expression for critical vibration amplitude needed to stabilize the drop a cr = gρ 1 ρ 2 (η 2 − η 1 ) [(ρ 1 + ρ 2 ) 2 + 2 (η 1 + η 2 ) 2 K −2 ω −2 ] ω 2 K(ρ 1 η 2 + 2ρ 2 η 1 )(ρ 2 1 − ρ 2 2 ) (5. 3) The stability boundary calculated during numerical simulation and the stability criterion (5.3) are plotted in Figure 7. One can see that for high frequency vibrations results of direct numerical simulations are in good agreement with the equation (5.3). For the case of low frequency our criterion slightly underestimate the critical vibration amplitude.
It is worth noting that for the case of stable drop results of calculation do not depend on calculation mesh. It is evident, that the settling rate increases in the case of stable inclusions. The similar results were obtained for rising drop which density is less than the density of outer liquid.

Conclusion
Numerical modeling of the liquid drop dynamics in porous medium saturated by the other liquid under gravity has been carried out. Simulations showed that the droplet is unstable and regardless of the ratio of liquid viscosities, the finger instability develops at the forefront of moving drop. Vibration effect on the plane displacement front in a porous medium has been analyzed in the framework of linear stability analysis and it is shown that the vibrations stabilizes the plane displacement front with respect to perturbation of short enough wavelength. Numerical modeling of vibration effect on the drop dynamics has demonstrated that the drop can be completely stabilized by vibrations of high intensity. The results of numerical modeling are in a good agreement with the obtained analytical criterion of the drop stability.