Towards Smoothed Particle Hydrodynamics Simulation of Viscous Fingering in Porous Media

Viscous fingering, which is observed at the interface between two immiscible fluids of different viscosities flowing through a porous medium when the more viscous fluid is displaced by the less viscous fluid [1,2], is a classical problem of fluid mechanics with important applications in oil recovery and earth drilling [36] and underpins the study of the wide range of Laplacian growth phenomena [7,8]. Although this phenomenon had been known to oil engineers for a long time, the first systematic studies of the interfacial instability during viscous fluid displacement appeared in the 1950s, using the Hele-Shaw cell as reference geometry based on the equivalence between the flow in a porous medium and the creeping flow between two parallel plates separated by a small gap [1,2]. Since then, viscous fingering is often referred to as SaffmanTaylor instability, which occurs when the displacement velocity exceeds a critical value:


Introduction
Viscous fingering, which is observed at the interface between two immiscible fluids of different viscosities flowing through a porous medium when the more viscous fluid is displaced by the less viscous fluid [1,2], is a classical problem of fluid mechanics with important applications in oil recovery and earth drilling [3][4][5][6] and underpins the study of the wide range of Laplacian growth phenomena [7,8]. Although this phenomenon had been known to oil engineers for a long time, the first systematic studies of the interfacial instability during viscous fluid displacement appeared in the 1950s, using the Hele-Shaw cell as reference geometry based on the equivalence between the flow in a porous medium and the creeping flow between two parallel plates separated by a small gap [1,2]. Since then, viscous fingering is often referred to as Saffman-Taylor instability, which occurs when the displacement velocity exceeds a critical value: Where 2 1 µ µ < . The problem of viscous fingering has been extensively studied form the experimental [9][10][11], theoretical [12][13][14], and computational point of view [15][16][17], for both Newtonian and non-Newtonian [18][19][20] fluids. Similar to other hydrodynamic instabilities, simulations of viscous fingering with commercial CFD codes using Eulerian meshes are difficult because it involves a moving interface usually characterized by large deformations. In this paper, Smoothed Particle Hydrodynamics (SPH) is used to simulate viscous fingering of Newtonian fluids in Hele-Shaw geometry. SPH is a fully Lagrangian computational technique, originally developed to simulate non-axisymmetric phenomena in astrophysics [21][22][23], that approximates the continuum fluid medium through the use of virtual particles and does not require a grid to calculate spatial derivate. In particular, a generic function at a given position r is represented as [21][22][23]: is the smoothing kernel and h is is the smoothing length that defines its influence region. The smoothing kernel can be an arbitrary function that (i) satisfies the normalisation condition, (ii) is identically zero outside the effective region defined by h, and (iii) tends to the Dirac delta when h→0. The integral representation of the spatial derivative of an arbitrary function reduces to the following equation: Introducing a smoothing function or smoothing kernel, the values of functions and spatial derivatives for a specific particle are approximated considering the interaction of that particle with a certain amount of neighbouring particles. This means that the physical quantities of a specific particle can be obtained by Through this approximation, the governing equations of fluid dynamics, i.e. the Navier-Stokes equations, are reduced to a set of ordinary differential equations with respect to time. Because of its Lagrangian nature, the SPH method has clear advantages over traditional grid base Eulerian methods for some fluid flow calculations, such as complex boundaries flows, multiphase fluids flows, free surfaces flows and non-Newtonian flows. In fact, since the particles simply follow their trajectories, fluid advection can be accomplished without difficulty. In particular, the absence of a fixed or adaptive mesh makes SPH particularly advantageous to track free surfaces and interfaces in comparison with, for example, with the Volume of Fluid method (VOF), where the exact position of interfaces is not determined a priori, and does not necessarily match the grid. However, the SPH method can be more computationally expensive than alternative techniques for a given application [24].
Historically, SPH took a relatively long time (decades) to become of widespread use in engineering fluid mechanics applications, mainly because in its original formulation some turbulent invariants are not conserved, and because sometimes the implementation of boundary conditions is not trivial; thus, several years were necessary to address these issues satisfactorily. Recently, the SPH method has been applied to simulate flow instabilities at the interface between two fluids, such as the Rayleigh-Taylor [25,26] and the Kelvin-Helmholtz instabilities [27,28], as well as the flow in porous media [29,30].However, the SPH approach has not been used to simulate the Saffman-Taylor instability so far.

Two-Phase SPH Algorithm
The SPH formulation of fluid dynamic equations is extensively discussed in the literature [21][22][23]31,32], and essentially consists in combining the conservation equations for mass, momentum and energy for the smoothed particles with a suitable constitutive equations that relates pressure with density. In most circumstances, a quasi-incompressible equation of state is used, so that the energy equation is not necessary. The particle density can be calculated either using a summation method: is the smoothing kernel of particle i evaluated at particle j (Eq. 10), or alternatively through the mass continuity equation: Where ij i j ν ν ν = − ; whilst Eq. (6) is more efficient than Eq (5) from the computational point of view, it does not ensure the mass conservation exactly [24]. The momentum equation can be written in the form of Newton's second law, and expresses the total force acting on a particle as the sum of a pressure force [23], a viscous force [24,33], and a body force (e.g., gravity): The choice of the constitutive equation and of the smoothing kernel can be somewhat controversial and depends on the characteristics of the flow under consideration [23,24,34]. Here, a quasi-incompressible flow equation of state [24] was chosen: Where c 2 is an artificial speed of sound, calculated as: In Eq. (9), v 0 and L 0 are the velocity and the length scales, respectively, and d is the density variation defined as Δρ⁄ρ.
A quantic spline kernel [31] was chosen as smoothing function: Where σ is a normalization constant that depends on the number of dimensions ( σ = 7/478 in 2D and σ = 3/359p in 3D, respectively),  the number of dimensions, R is the ratio between the modulus of the distance vector, r and the smoothing length h. Although this choice is computationally more expensive than a cubic spline kernel, it is more stable for flows with low Reynolds numbers [24]. Finally, the time step was selected small enough not only to satisfy the CFL condition, but also to resolve adequately the particle acceleration and the viscous diffusion: The main idea behind the SPH algorithm is to solve the Poiseuille problem between two infinite parallel plates where two different fluids (represented by two distinct sets of particles) are initially placed side by side, as shown in Figure 1. This means that, instead of calculating only the interactions between pairs of fluid particles and those between fluid particles and boundary particles, the algorithm must evaluate the interactions between pairs of fluid particles of the same set, the interactions between pairs of fluid particles belonging to different sets, and the interactions between particles of each set and boundary particles.
In standard SPH algorithms, the properties of boundary particles are calculated based on the nearest fluid particle; however, this can generate significant errors when boundary particles interact with fluid particles belonging to different sets. To overcome this issue, two overlapping sets of boundary particles were used, one for each set of fluid particles.

Validation
The code was initially validated against the analytical solutions of the single-phase, time-dependent Poiseuille flow; the validation was then repeated placing in the computational domain two fluids initially separated, identified by different sets of particles but with the same density and viscosity. A viscous fluid initially at rest starts to flow between two parallel plates located at y=0 and y=L, driven by a body force F(i.e. the pressure gradient) parallel to the plates (x-axis), until it reaches the well-known steady-state solution for the velocity profile: Where F is the body force per unit mass and  the kinematic viscosity. The full transient solution is [24]: In this simulation, the following parameters were chosen: F in the x-direction, while the two plates were modelled with 105 particles each, 5 in the y-direction and 21 in the x-direction. The initial smoothing length was taken equal to 1.33 times the initial particle space, which in turn was calculated as one twentieth of the distance between the two plates; the time step was set to 0.0001s. The simulation reached a steady state after around 6000 time steps however, in order to check the actual effect of the periodic boundary, a plot of the particle distribution was taken after 50,000 time steps and displayed in Figure 2. The comparison of the results SPH simulation with the analytical solution (Eq. 13), displayed in Figure 2, shows a very close agreement, with an average relative error of 1.4% after 100 time steps and 0.3% after 1000 and 6000 time steps. The code was then tested using two fluids characterised by the same density and viscosity. In this way, the velocity profile obtained must be identical to the velocity profile of the single-phase Poiseuille flow. All parameters, including the initial smoothing length and the time step length, were the same used previously. The problem domain was a rectangle of dimensions 0.002m x 0.001m, with the first fluid placed in the half on left side. The latter was modelled with 380 particles, 19 in the y-direction and 20 in the x-direction while the second fluid was modelled with 399 particles, 19 in the y-direction and 21 in the x-direction. The first boundary was modelled with 200 particles, 10 in the y-direction and 20 in the x-direction and the second boundary was modelled 210 particles, 10 in the y-direction and 21 in the x-direction. The same body force was applied to both fluids particles, and the resulting velocity profile obtained was exactly the same as the one shown in the Figure 2b, while the particle distribution after 50000 time steps is shown in the Figure 3.

Viscous Fingering Simulations
SPH simulations of viscous fingering were conducted using two fictitious fluids characterised by the same value of density and different values of the kinematic viscosity. The channel containing the two fluids was a rectangle of length 5x10 -3 m and width 10 -3 m, with the lower-viscosity fluid ("fluid_1") placed in the left half of the channel and pushing the higher-viscosity fluid ("fluid_2") initially at rest in the right half of the channel, as shown in Figure 1. Fluid_1 was modelled with 950 particles, 19 in the y-direction and 50 in the x-direction, while fluid_2 was modelled with 969 particles, 19 in the y-direction and 51 in the x-direction. Boundaries were modelled with 500 particles for fluid_1, 10 in the y-direction and 50 in the x-direction, and with 510 particles, 10 in the y-direction and 51 in the x-direction, for fluid_2.
The reference density of the two fluids was set to 1000kg/m 3 , the kinematic viscosity of fluid_1 was kept constant at a value of 10 -6 m 2 /s while the value of the viscosity of fluid_2 was 10 -5 m 2 /s. The Reynolds number was varied in the range between 0.1 and 1; this condition was implemented by calculating the laminar pressure gradient corresponding to the Reynolds number value, and applying it as a body force acting on particles corresponding to both fluids. The upper limit of the Reynolds number magnitude was determined by the CFL condition (Eq. 11), using a time step of 10 -4 s for all simulations. The initial smoothing length was equal to 1.33 times the initial particles spacing, which in turn was calculated as one twentieth of the channel width. Simulations were run for 10,000 time steps, corresponding to 10s.

Conclusion
A two-phase mesh-free Lagrangian SPH code was developed and validated to simulate the displacement of fluids in porous media by means of a fluid having different viscosity. The code features an original implementation of boundary particles, which create a virtual channel for each fluid to avoid discontinuities at the contact point where the interface between fluids meets the channel wall. Preliminary results show that SPH can reproduce the formation of a low-viscosity finger penetrating into the higher-viscosity fluid during displacement. Because of its intrinsically Lagrangian, meshfree nature, SPH is a promising method to simulate viscous fingering in complex geometries; in addition, it can easily incorporate non-Newtonian constitutive equations to account for shear-thinning and/or viscoplastic behaviours.