Crimson Publishers Publish With Us Reprints e-Books Video articles

Full Text

Progress in Petrochemical Science

Towards Smoothed Particle Hydrodynamics Simulation of Viscous Fingering in Porous Media

Ronchi N and Bertola V*

Laboratory of Technical Physics, University of Liverpool, UK

*Corresponding author: Bertola V, Laboratory of Technical Physics, School of Engineering, University of Liverpool, Brownlow Hill, UK

Submission: May 05, 2018;Published: May 10, 2018

DOI: 10.31031/PPS.2018.01.000523

ISSN 2637-8035
Volume1 Issue5


The mesh-free Lagrangian Smoothed Particle Hydrodynamics (SPH) method is used to simulate the problem of viscous fingering in Hele-Shaw geometry. Viscous fingering occurs when a lower-viscosity fluid displaces a higher-viscosity fluid in a narrow channel, and is a major concern in the removal of drilling mud’s from oil well bores and in secondary and tertiary oil recovery. The flow field was modelled using two sets of particles with different properties, initially placed in the left half and in the right half of the channel, respectively; in particular, the two fluids had the same density and a viscosity ratio of 1:10. 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, mesh-free 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.


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- 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-11], theoretical [12- 14], and computational point of view [15-17], for both Newtonian and non-Newtonian [18-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-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-23]:

Where W(r − r′, h) 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 summing the relevant properties of all the particles which lie within the range of the smoothing function, whose values are smoothed by the kernel function itself:

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-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:

Where Wij =W ri − rj h is the smoothing kernel of particle i evaluated at particle j (Eq. 10), or alternatively through the mass continuity equation:

Where νiji −ν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 c2 is an artificial speed of sound, calculated as:

In Eq. (9), v0 and L0 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.

figure 1:Initial particle distribution of the two fluids and boundary particles.


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 = 2.10-4 m/s2, L = 1.10-3m, v = 1.10-6m2/s so that the Reynolds number based on the fluid maximum velocity was 2.5•10-2. The problem domain was a square of 0.001m x 0.001m. The fluid was modelled with 399 particles, 19 in the y-direction and 21 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.

figure 2:Steady-state (t=5s) particle distribution in single-phase Poiseuille flow (A), and comparison of computed transient velocity profiles with the analytical solution (Eq. 13) (B).

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.

figure 3:Steady-state particle distribution obtained for the Poiseuille flow of two fluids with identical properties.

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-3m and width 10-3m, 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/m3, the kinematic viscosity of fluid_1 was kept constant at a value of 10-6m2/s while the value of the viscosity of fluid_2 was 10-5m2/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.

figure 4:Evolution of viscous fingering in fluids with kinematic viscosity ratio 2/1=10 and Reynolds number Re=0.1, at t=0.5s (A), t=5s (B), t=10s (C).

Figure 4-6 display the evolution of particle distribution at three different Reynolds numbers (Re=0.1, Figure 4; Re=0.5, Figure 5; Re=1, Figure 6). Unlike in the case of fluids with identical properties, one can observe an evolution of the interface that indicates the development of viscous fingering. Remarkably, fingering is obtained purely as a result of interactions among particles, without explicit modelling the interfacial tension.

figure 5:Evolution of viscous fingering in fluids with kinematic viscosity ratio 2/11=10 and Reynolds number Re=0.5, at t=0.5s (A), t=5s (B), t=10s (C).

figure 6:Evolution of viscous fingering in fluids with kinematic viscosity ratio 2/11=10 and Reynolds number Re=1, at t=0.5s (A), t=5s (B), t=10s (C).


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.


  1. Lewis DJ (1950) The instability of liquid surfaces when accelerated in a direction perpendicular to their planes. II. Proc R Soc Lond 202(1068): 81-96.
  2. Saffman PG, Taylor G (1958) The penetration of a fluid into a porous medium or Hele-Shaw cell containing a more viscous liquid. Proc R Soc Lond 245(1242): 312-329.
  3. Latil M (1980) Enhanced oil recovery. Editions OPHRYS.
  4. Thomas S (2008) Enhanced oil recovery-an overview. Oil & Gas Science and Technology-Rev. IFP 63(1): 9-19.
  5. Bittleston S, Ferguson J, Frigaard I (2002) Mud removal and cement placement during primary cementing of an oil well-Laminar non-Newtonian displacements in an eccentric annular Hele-Shaw cell. Journal of Engineering Mathematics 43(2-4): 229-253.
  6. Pelipenko S, Frigaard I (2004) On steady state displacements in primary cementing of an oil well. J Eng Math 46(1): 1-26.
  7. Gustafsson B, Vasilev A (2006) Conformal and potential analysis in Hele- Shaw cells. Springer.
  8. Vasilev A (2009) From the Hele-Shaw Experiment to Integrable Systems: A Historical Overview. Complex Analysis and Operator Theory 3(2): 551-585.
  9. Pitts E (1980) Penetration of fluid into a Hele-Shaw cell: The Saffman- Taylor experiment. Journal of fluid mechanics 97(1): 53-64.
  10. Couder Y, Gerard N, Rabaud M (1986) Narrow fingers in the Saffman-Taylor instability. Phys Rev A Gen Phys 34(6): 5175-5178.
  11. Tabeling P, Zocchi G, Libchaber A (1987) An experimental study of the Saffman-Taylor instability. Journal of Fluid Mechanics 177: 67-82.
  12. David AK, Levine H (1986) Theory of the Saffman-Taylor ‘‘finger’’ pattern I. Phys Rev A 33: 2621.
  13. David AK, Levine H (1986) Theory of the Saffman-Taylor ‘‘finger’’ pattern II. Phys Rev A Gen Phys 33: 2634.
  14. Hakim CR, Dombre V, Pomeau T Y, Pumir A (1988) Analytic theory of the Saffman-Taylor fingers. Phys Rev A Gen Phys 37(4): 1270-1283.
  15. Gretar T, Hassan A (1983) Numerical experiments on Hele Shaw flow with a sharp interface. Journal of Fluid Mechanics 136: 1-30.
  16. Degregoria AJ, Schwartz LW (1986) A boundary-integral method for two-phase displacement in Hele-Shaw cells. Journal of Fluid Mechanics 164: 383-400.
  17. Whitaker N (1997) Some numerical methods for the hele-shaw equations. Journal of Computational Physics 111(1): 81-88.
  18. Wilson SDR (1990) The Taylor-Saffman problem for a non-Newtonian liquid. Journal of Fluid Mechanics 220: 413-425.
  19. Mora S, Manna M (2009) Saffman-Taylor instability for generalized newtonian fluids. Phys Rev E Stat Nonlin Soft Matter Phys 80(1 Pt 2): 016308.
  20. Mora S, Manna M (2010) Saffman-Taylor instability of viscoelastic fluids: From viscous fingering to elastic fractures. Phys. Rev E 81(2).
  21. Lucy LB (1977) A numerical approach to the testing of the fission hypothesis. The astronomical journal 82: 1013-1024.
  22. Gingold RA, Monaghan JJ (1977) Smoothed particle hydrodynamics-theory and application to non-spherical stars. Monthly notices of the royal astronomical society 181: 375-389.
  23. Monaghan JJ (1992) Smoothed particle hydrodynamics. Annual review of astronomy and astrophysics 30: 543-574.
  24. Morris JP, Fox PJ, Zhu Y (1997) Modeling low Reynolds number incompressible flows using SPH. Journal of computational physics 136(1): 214-226.
  25. Rahmat A, Tofighi N, Shadloo M, Yildiz M (2014) Numerical simulation of wall bounded and electrically excited Rayleigh–Taylor instability using incompressible smoothed particle hydrodynamics. Colloids and Surfaces A: Physicochemical and Engineering Aspects 460: 60-70.
  26. Tartakovsky AM, Meakin P (2005) A smoothed particle hydrodynamics model for miscible flow in three-dimensional fractures and the two-dimensional Rayleigh-Taylor instability. Journal of Computational Physics 207(2): 610-624.
  27. Price DJ (2008) Modelling discontinuities and Kelvin–Helmholtz instabilities in SPH. Journal of Computational Physics 227(24): 10040-10057.
  28. Shadloo MS, Yildiz, M (2011) Numerical modeling of Kelvin-Helmholtz instability using smoothed particle hydrodynamics. International Journal for Numerical Methods in Engineering 87(10): 988-1006.
  29. Holmes DW, Williams JR, Tilke P (2011) Smooth particle hydrodynamics simulations of low Reynolds number flows through porous media. International Journal for Numerical and Analytical Methods in Geomechanics 35(4): 419-437.
  30. Tartakovsky AM, Meakin P (2006) Pore scale modeling of immiscible and miscible fluid flows using smoothed particle hydrodynamics. Advances in water resources 29(10): 1464-1478.
  31. Liu GR, Liu M (2003) Smoothed particle hydrodynamics: a meshfree particle method. World Scientific.
  32. Liu M, Liu G (2010) Smoothed particle hydrodynamics (SPH): an overview and recent developments. Archives of computational methods in engineering 17(1): 25-76.
  33. Monaghan J (1995) Heat conduction with discontinuous conductivity. Applied Mathematics Reports and Preprints 95(18): 7.
  34. Monaghan J, Rafiee A (2013) A simple SPH algorithm for multi‐fluid flow with high density ratios. International Journal for Numerical Methods in Fluids 71(5): 537-561.

© 2018 Bertola V. This is an open access article distributed under the terms of the Creative Commons Attribution License , which permits unrestricted use, distribution, and build upon your work non-commercially.