Study of Multi-phase Flow in Porous Media: Comparison of SPH Simulations with Micro-model Experiments

Abstract

We present simulations and experiments of drainage processes in a micro-model. A direct numerical simulation is introduced which is capable of describing wetting phenomena on the pore scale. A numerical smoothed particle hydrodynamics model was developed and used to simulate the two-phase flow of immiscible fluids. The experiments were performed in a micro-model which allows the visualization of interface propagation in detail. We compare the experiments and simulations of a quasistatic drainage process and pure dynamic drainage processes. For both, simulation and experiment, the interfacial area and the pressure at the inflow and outflow are tracked. The capillary pressure during the dynamic drainage process was determined by image analysis.

1 Introduction

In this work, we present a comparison between direct numerical simulations of immiscible two-phase flow in a porous structure based on smoothed particle hydrodynamics (SPH) and the results obtained from micro-model experiments.

First SPH simulations of multi-phase flow in porous structures on, the pore scale have been performed by Tartakovsky and Meakin (2006) and Tartakovsky et al. (2015). Still comparisons with experiments are very rare and focus mostly on the dynamics of a single interface, like a capillary rise. The need for an experimental work where a number of capillaries, like a porous structure, would be used became evident.

Both quasi-static and dynamic drainage experiments were performed in a poly-dimethyl-siloxane (PDMS) micro-model (Karadimitriou et al. 2012). The fluids used as the wetting and the non-wetting phase were Fluorinert FC-43 and water dyed with ink, respectively. The dynamic drainage experiments were performed under constant pressure conditions. Flow through the micro-model was visualized with the use of a custom open-air microscope. The acquired images were stored in a computer for further processing (Karadimitriou et al. 2013).

The images acquired by the visualization setup were subsequently processed to obtain information on the interfacial area between the wetting and the non-wetting phase, as well as between the two fluids and the solid phase. Phase saturation and the pressure in the reservoir at the inlet were also captured. In this way, a direct comparison of experiments and simulations was obtained.

For the simulations a modified version of the Lagrangian SPH method (Gingold and Monaghan 1977) is used to easily track the moving interfaces of the multi-phase flow. To account for interface and contact line phenomena, the Navier–Stokes equations are extended by the Continuum Surface Force model (CSF) (Brackbill et al. 1992) to describe the forces acting on the interface and the contact line force model (CLF) (Huber et al. 2013) to describe wall–fluid–fluid interactions. These models allows us to track the extension of all interfaces and contact lines.

The simulations are performed in 2D to capture the details of the horizontal interface propagation. The influence of the top and bottom wall of the micro-model on the viscous drag is included into the Navier–Stokes equation. Additionally, a constant contribution to the capillary pressure is assumed by the interfacial curvature in vertical direction.

Fig. 1

Micro-model pore network

2 Experiment

2.1 Experimental Setup

The experimental setup is the one used in various studies by Karadimitriou et al. (2012, 2013). The micro-model used for the experiment had a rectangular shape with dimensions $2.5\times 1\,\hbox {mm}$. It consists of a big void, which corresponds to the size of the micro-model, with a random distribution of round pillars. This configuration creates a pore network with average pore diameter of $160\,\upmu \hbox {m}$ (Fig. 1). The micro-model was fabricated according to the work of Karadimitriou et al. (2013). First a silicon wafer with the imprinted micro-model was created using photolithography. Subsequently, this wafer was used as mold on one PDMS slab. A second PDMS slab was then used to cover the micro-model.

The experimental setup consists of three major parts. The first part is the image acquisition system, which consists of a light source, a prism, camera lenses, and a monochromatic camera. The second part, consists of two reservoirs containing the two liquid phases where an equal number of pressure transducers are fixed to measure the pressure difference in the reservoirs (Fig. 2). The inlet and outlet of the micro-model are connected to the two reservoirs. The third part of the setup contains the pumping system, which controls the applied pressure to the reservoir of the non-wetting phase. More details about the experimental setup can be found in the work of Karadimitriou et al. (2013).

Fig. 2

Micro-model position on the setup

2.2 Two-Phase Flow Experiments

The two-phase flow experiment has been performed according to Karadimitriou et al. (2013). The micro-model is connected with the two reservoirs. Each reservoir contains a different liquid phase, one being the wetting phase and one the non-wetting phase. In this experiment, due to the hydrophobicity of PDMS, water (blue ink) is the non-wetting phase and Fluorinert FC-43 is the wetting phase. Figure 1 shows the micro-model fully wetted with Fluorinert. In all experiments, we acquired images of the phase propagation. From those images, all interfacial areas and the capillary pressure at the wn interface was determined.

2.2.1 Dynamic Drainage Experiment

The dynamic drainage experiment consisted of one full drainage process. During drainage a constant external pressure was applied to the reservoir containing the non-wetting phase. The pressure is controlled by a pressure valve (Bronkhorst) and maintained a stable pressure during the experiment. In order to achieve dynamic conditions, the applied pressure was kept at $1.86\,\text {kPa}$ during the drainage process. Images of the process were taken with a rate of 5 fps.

2.2.2 Quasistatic Drainage Experiment

The quasistatic experiment has been performed using a flow control. The drainage process was paused for eight times to let the system relax. One image was taken in each relaxed condition, and the current pressure in the reservoir was measured.

3 Model Description

Isothermal flow of immiscible and incompressible Newtonian fluids in a continuous domain is described by the Navier–Stokes equations. The Lagrangian formulation reads

$$\begin{aligned} \rho \,\frac{\mathrm{D}\mathbf {v}}{\mathrm{D}t}=-\nabla p + \mu \nabla ^2\mathbf {v} + \rho \,\mathbf {g} + \mathbf {F}_{wn}^\mathrm{VOL} + \mathbf {F}_{wns}^\mathrm{VOL}, \end{aligned}$$

(1)

where the term $-\nabla p + \mu \nabla ^2\mathbf {v}$ results from a pressure gradient and viscous stress in the system. Further $\mathbf {g}$ represents any body force like gravity. $\mathbf {F}_{wn}^\mathrm{VOL}$ and $\mathbf {F}_{wns}^\mathrm{VOL}$ are volumetric forces which represent the forces acting at the interface between the wetting and non-wetting phase and the forces at the contact line between the two fluids and the wall. For incompressible fluids, the continuity equation simplifies to

$$\begin{aligned} \frac{\mathrm{D}\rho }{\mathrm{D}t}=-\rho \,(\nabla \cdot \mathbf {v}) = 0. \end{aligned}$$

(2)

3.1 Implementation in SPH

We chose the SPH method because of its mesh-free nature. Hence SPH can treat the propagation of multiple interfaces in a pure Lagrangian manner. Also the splitting and merging of existing interfaces is an easy task in SPH. A comparable accuracy of SPH with other methods like lattice Boltzmann or FEM in flow through porous media was shown by Sivanesapillai et al. (2014). SPH is an interpolation method where properties are evaluated at certain interpolation points to approximate a continuous field. These points, which represent a certain volume, are called particles and move in the system in a Lagrangian manner. In SPH the interpolation formula for any quantity $A(\mathbf {r})$ is an approximation to an integral interpolant of the form (Monaghan 2011)

$$\begin{aligned} A(\mathbf {r}) = \int A(\mathbf {r}')\,W(\mathbf {r}-\mathbf {r}',h)\,\mathbf {dr}' \approx \sum _j \frac{m_j\,A_j}{\rho _j}W_{ij}, \end{aligned}$$

(3)

where $A_j = A(\mathbf {r}_j)$ and $W_{ij} = W(\mathbf {r}_i-\mathbf {r}_j,h)$ is a kernel function with smoothing length $h=2.1$. We use the C2 spline function by Wendland (1995)

$$\begin{aligned} W(\mathbf {r},h) = \frac{7}{4\pi h^2}{\left\{ \begin{array}{ll} (1-\frac{q}{2})^4\,(2q+1) &{}\text {if}\,\, q<2\\ 0 &{}\text {else} \end{array}\right. }, \end{aligned}$$

(4)

where $q = |\mathbf {r}|/h$. This kernel was shown to have a lower tendency for particle clustering compared to other kernels (Szewc et al. 2012). A detailed introduction into the SPH theory can be found for example in the review of Monaghan (2005).

3.1.1 SPH Formulation for Multi-phase Flow

Our multi-phase model is based on the model presented by Szewc et al. (2012), with the exception of the viscous force, where we chose the formulation of Morris et al. (1997). We chose this formulation, because it produces an accurate velocity profile in low Reynolds number flow in a channel (Morris et al. 1997). The acceleration due to the viscous force is formulated as

$$\begin{aligned} \left( \left( \frac{1}{\rho }\nabla \cdot \mu \nabla \right) \mathbf {v}\right) _i=\sum _j \frac{m_j\,(\mu _i+\mu _j)\,\mathbf {r}_{ij}\cdot \nabla _i W_{ij}(h)}{\rho _i\,\rho _j\mathbf {r}_{ij}^2}\mathbf {v}_{ij}, \end{aligned}$$

(5)

with $\mathbf {r}_{ij} = \mathbf {r}_{i} - \mathbf {r}_{j}$ .This formulation conserves linear momentum exactly, while angular momentum is only approximately conserved. The density of the particles is assumed to be constant $\rho _i = m_i / V_0$ because of incompressibility. Here $V_0$ is the volume of one particle, which is constant for all the particles. The pressure force is originally introduced by Morris (2000) as

$$\begin{aligned} -\left( \frac{1}{\rho } \nabla p\right) _i = -\sum _j m_j \frac{p_i + p_j}{\rho _i \, \rho _j} \nabla _i W\left( \mathbf {r}_{ij},h\right) . \end{aligned}$$

(6)

3.1.2 Surface Tension in SPH

Different numerical methods are available to model surface tension in SPH. The most common one is the CSF model which was introduced by Brackbill et al. (1992) and brought to SPH by Morris (2000). This approach does not conserve momentum exactly and some spurious currents occur at the interface. However, the method has proven to be very stable (Szewc et al. 2013). An extension to the CSF model which includes forces at the contact line (CLF) was proposed by Huber et al. (2013). There is a momentum conserving variant of the CSF method available, where the effect of surface tension is added as a correction to the momentum stress tensor (Lafaurie et al. 1994). We tested the implementation of Hu and Adams (2006), but the effect of spurious currents, especially at the contact line, was larger compared to the classical CSF model. Tartakovsky and Meakin (2006) used a less computational demanding method based on a methodically different approach from Nugent and Posch (2000). They describe the macroscopic effect of surface tension with attractive and repulsive forces to avoid the calculation of the local curvature at the interface. One disadvantage of this method is that the force will not vanish at straight interfaces. Considering all methods we decided to use the CSF/CLF model for this work. With the assumption of a constant surface tension $\sigma _{wn}$ between wetting and non-wetting phase, the CSF model reads

$$\begin{aligned} \mathbf {F}_{wn}^\mathrm{VOL} = \mathbf {f}_{wn}\,\delta _{wn} = \sigma _{wn}\,\kappa _{wn}\,\mathbf {n}_{wn}, \end{aligned}$$

(7)

where $\mathbf {f}_{wn}$, $\delta _{wn}$, $\sigma _{wn}$, $\kappa _{wn}$ and $\mathbf {n}_{wn}$ are the surface force, volume reformulation parameter, surface tension coefficient, curvature and the normal vector of the interface between the two immiscible phases. The volume reformulation is necessary to transfer the forces acting at the interface or contact line into the Navier–Stokes equations for the fluid domain. The normal vector is calculated using the color (or indicator) function c.

$$\begin{aligned} \mathbf {n}_{wn} = \frac{\nabla c}{[c]} \end{aligned}$$

(8)

where $[c]= |c_w-c_n|$ is the change in c at the interface. In SPH Eq. (8) can be written as

$$\begin{aligned} \mathbf {n}_{wn,i}&= \frac{1}{[c]}\sum _j\frac{m_j}{\rho _j}\left( c_j-c_i \right) \nabla _{i}W_{ij}, \end{aligned}$$

(9)

Here $|\mathbf {n}_{wn}|$ can be used as volume reformulation $\delta _{wn,i}$. The local curvature of the interface is defined as

$$\begin{aligned} \kappa _{wn} = -\nabla \cdot \hat{\mathbf {n}}_{wn}, \end{aligned}$$

(10)

where $\hat{\mathbf {n}}_{wn}$ is the normalized normal vector. The curvature is discretized in SPH as:

$$\begin{aligned} \kappa _{wn,i} = -\sum _j \frac{m_j}{\rho _j}\left( \hat{\mathbf {n}}_{wn,j} - \hat{\mathbf {n}}_{wn,i}\right) \cdot \nabla _i W_{ij} \end{aligned}$$

(11)

We also tested other implementations for the CSF model. Adami et al. (2010b) proposed a formulation which performs the volume reformulation with a density weighting. This formulation tends to an increased mixing of the two phases. We also implemented the formulation of Hu and Adams (2006) which results in increased spurious currents at the contact line. The chosen formulation still showed a tendency for mixing of the two phases and spurious currents at the interphase. This is a common issue for the CSF model with a high ratio between surface tension and viscosity. However, the mixing effect was small enough to avoid an extra repulsive force at the interface, as the one proposed by Grenier et al. (2009).

The momentum balance at a contact line tangential to the wall including the volume reformulation reads (Huber et al. 2013; Hassanizadeh and Gray 1993)

$$\begin{aligned} \mathbf {F}_{wns}^\mathrm{VOL} \cdot \varvec{\nu }_{ns} = \bigg (\sigma _{ns} - \sigma _{ws} + \sigma _{wn}\,(\underbrace{\hat{\varvec{\nu }}_{ns} \cdot \hat{\varvec{\nu }}_{wn}}_{-\cos {\theta }})\bigg )\,\delta _{wns}, \end{aligned}$$

(12)

where $\hat{\varvec{\nu }}_{kl}$ are unit vectors pointing away from the contact line (see Fig. 3a). They lie on the kl-interface and are orthogonal to the contact line itself. $\delta _{wns}$ is the volume reformulation parameter for the CLF (Huber et al. 2013).

3.2 Corrected SPH

We use the correction of the kernel and the correction of the gradient introduced by Bonet and Lok (1999) in this research work for all computed properties. We do this for three reasons:

• 1. The calculation of the normal vector $\mathbf {n}$ and the curvature $\kappa$ does not have full support in the domain (Morris 2000; Adami et al. 2010a) (see Fig. 3b, c).

• 2. The simulations are performed with open boundaries. This means that new particles can enter the system and others leave it. In this inlet and outlet zone, the particle distribution is not perfect. In order to minimize the errors at the open boundaries the application of corrected kernels and their gradients is crucial.

• 3. The convergence with increasing particle resolution in a disturbed system is only guaranteed with higher-order methods (Fatehi and Manzari 2011; Tartakovsky et al. 2015).

Fig. 3

a The contact line at the connection point of the unit vectors $\hat{\varvec{\nu }}_{kl}$, where each vector $\hat{\varvec{\nu }}_{kl}$ is located tangential on the corresponding interface. b The missing support in the calculation of $\mathbf {n}_i$ close to a wall and in c the curvature $\kappa _i$ is only calculated using particles which hold the condition $|\mathbf {n}_j|> 0$

The corrected kernel can be described by:

$$\begin{aligned} \tilde{W}_{ij} = \frac{W_{ij}}{\sum _{j}V_jW_{ij}} = \frac{W_{ij}}{S_i}, \end{aligned}$$

(13)

where $S_i$ is called the Shepard function of particle i. The gradient of this corrected kernel becomes

$$\begin{aligned} \nabla \tilde{W}_{ij} = \nabla \left( \frac{W_{ij}}{S_i}\right) = \frac{\nabla W_{ij} - \frac{W_{ij}\,\nabla S_{i}}{S_{i}}}{{S_{i}}}, \end{aligned}$$

(14)

with $\nabla S_{i} = {\sum }_j V_j \nabla W_{ij}$. To ensure that any linear property field is evaluated correctly and angular momentum is conserved, the following condition must be fulfilled:

$$\begin{aligned} \mathbf {L}_i = \bigg (\sum _j V_j \nabla \tilde{W}_{ij} \otimes (-\mathbf {r}_{ij}) \bigg )^{-1}. \end{aligned}$$

(15)

Thus the corrected gradient of the corrected kernel reads:

$$\begin{aligned} \tilde{\nabla } \tilde{W}_{ij} = \mathbf {L}_i \nabla \tilde{W}_{ij} \end{aligned}$$

(16)

Those corrected gradients of the corrected kernels were used for all formulations.

3.3 Incompressibility Treatment and Time Integration

In our work, we use the truly incompressible SPH (ISPH) method based on the projection method introduced by Cummins and Rudman (1999). In ISPH the pressure can be calculated by solving a pressure Poisson equation (PPE). There are two criteria available from the continuity equation (Eq. 2) to satisfy the incompressibility condition. Either the divergence of the velocity is zero or the density-field is invariant over one time step.

$$\begin{aligned} \nabla \cdot \left( \frac{1}{\rho } \nabla p\right) = \frac{\mathrm{div} \, {\mathbf {v}^{*}}}{\Delta t} = \frac{\rho _0-\rho ^*}{\rho _0 \Delta t^2}. \end{aligned}$$

(17)

Equation 17 is discretized as

$$\begin{aligned} \nabla \cdot \left( \frac{1}{\rho } \nabla p\right) _i&= \sum _j \frac{m_j}{\rho _j} \frac{4}{\rho _i + \rho _j} \frac{p_{ij} \mathbf {r_{ij}} }{|\mathbf {r_{ij}}|^2} \cdot \nabla _i W_{ij}\left( h\right) \end{aligned}$$

(18)

$$\begin{aligned} \left( \mathrm{div} \, {\mathbf {v}^{*}} \right) _i&= \sum _j \frac{m_j}{\rho _j}\left( \mathbf {v_{j}}^* - \mathbf {v_{i}}^*\right) \cdot \nabla _i W_{ij}\left( h\right) \end{aligned}$$

(19)

The intermediate velocity $\mathbf {v}^*$ is computed from the momentum balance, but excluding the pressure force:

$$\begin{aligned} \mathbf {v}^*=\left( \frac{\mu }{\rho }\nabla ^2\mathbf {v} + \mathbf {g} + \frac{1}{\rho }\mathbf {F}_{wn}^\mathrm{VOL} + \frac{1}{\rho }\mathbf {F}_{wns}^\mathrm{VOL}\right) \Delta t + \mathbf {v}^t. \end{aligned}$$

(20)

We use a time integration scheme with a divergence-free velocity field and a density-invariant field (ISPH_DFDI) introduced by Hu and Adams (2007). The velocity at the new time step is computed by:

$$\begin{aligned} \mathbf {v}^{t+1} = \mathbf {v}^* - \left( \frac{1}{\rho } \nabla p \right) \cdot \Delta t \end{aligned}$$

(21)

We additionally use the particle shifting technique (ISPH_DFS) from Xu et al. (2009). The shifting was implemented at the end of the time step to overcome tensile instabilities which could for instance occur from a homogeneous Cartesian initial particle distribution or the change in interfacial area. The integration scheme can be found in the work of Xu et al. (2009).

3.4 Boundary Conditions

3.4.1 Solid Wall Boundary Conditions

The imbibition and drainage processes in the micro-model happen at low Reynolds numbers. Additionally, the Knudsen numbers of the two present liquids in the smallest channels is still very small ($Kn<10^{-5}$). Therefore, it is valid to use no-slip boundary conditions introduced by Morris (2000) for both liquids in the SPH simulation to realistically model laminar flow. It is important to mention, that there is no projected velocity at the wall particles when the pressure Poisson equation is solved (Eq. 19). The no-slip condition is only applied when calculating the viscous force (see Fig. 4a).

Fig. 4

a The artificial velocity for boundary particles to simulate a no-slip boundary condition for a particle at a solid wall. b The velocities of the mirrored particles at the open boundary

3.4.2 Inflow/Outflow Condition

In this work, we use Dirichlet boundary conditions of the projected pressure field at the inlet and outlet. We impose a divergence-free boundary condition at the inlet and outlet by mirroring the particles at the open boundary and keeping the velocities of the mirrored particles. In Fig. 4b the mirrored particles at the open boundary are shown. Also the occurrence of small void spaces at the inlet and outlet is represented by the different distance of the particles from the boundary. The same mirroring technique is used to specify the Dirichlet boundary condition for the pressure. This condition at the domain boundary (where no particles are located) is fulfilled by defining the pressure of the mirrored particles $j'$. For a mirrored particle $j'$ this pressure is defined as:

$$\begin{aligned} p'_j = 2 p_w - p_j , \end{aligned}$$

(22)

where $p_w$ is the desired pressure at the inlet or outlet. The conditions for open pressure boundary conditions then are given as:

$$\begin{aligned} \begin{aligned} \mathrm{div} \, \mathbf {u} |_w = 0 \\ p_w = p_0, \end{aligned} \end{aligned}$$

(23)

3.5 Implementation

The presented SPH models are implemented in SiPER.¹ SiPER is an ISPH code developed by our group at the University of Stuttgart. It is specialized on multi-phase flow in complex porous media. It is a MPI parallel implementation in standard C language. We use the boomerAMG preconditioner and BiCGStab solver from HYPRE and PETSc library for solving the PPE.

4 Simulation

4.1 Validation of the 2D Simulation

The imbibition and drainage process in a porous structure is driven by interfacial forces and a pressure gradient, respectively. Moreover the propagation speed is constrained by viscous forces. Before performing the simulations of the micro-model, the model was validated by two test cases. In the first test case, the flow profile in a single-phase Poiseuille flow is validated and in the second test case a dynamic capillary rise is compared with a reduced model.

4.1.1 Poiseuille Flow

The Poiseuille flow is a simple test case to verify a correct implementation of no-slip boundaries at the solid wall and viscous forces in general. In most literature of SPH, Poiseuille flow is applied by using periodic boundary conditions in the direction of the flow and applying a body force to accelerate the fluid. In this work we apply a pressure gradient between inlet and outlet with Dirichlet boundary conditions for the pressure as driving force (see Sect. 3.4.2). In this way we can use the Poiseuille flow also as test case for the pressure field and the acceleration term arising from a pressure gradient.

We consider a two-dimensional, rectangular domain with a length of $L_x = 3D = 480\,\upmu \hbox {m}$ and a height of $L_y = D$. In y-direction we apply no-slip wall boundary conditions. The fluid has a density of $\rho = 1000\,\hbox {kg}/\hbox {m}^3$ and a dynamic viscosity of $\mu = 0.001\,\text {Pa s}$. The Reynolds number is $Re = \frac{\rho \mathbf {u}_\mathrm{avg} D}{\mu } = 7.1$. The dimensionless velocity is defined as $\mathbf {u}^* = \frac{\mathbf {u}}{\mathbf {u}_\mathrm{max}}$. The dimensionless positions are $x^* = \frac{x}{D}$ and $y^* = \frac{y}{D}$. The corresponding pressure difference between in- and outlet is $\Delta p^* = \frac{\Delta p}{\mathbf {u}_\mathrm{avg}^2 \rho } = 5.06$. Initially the system is at rest. The domain size is set to 7500 particles.

In Fig. 5a velocity profiles for three times $t^* = \frac{t \cdot \mathbf {u}_\mathrm{max}}{D}$ are shown. In Fig. 5b, the error in the pressure field is plotted over the dimensionless position $x^*$. It is computed by $\hbox {err} = \frac{p_\mathrm{simulation}-p_\mathrm{analytic}}{p_\mathrm{anayltic}}$. The new introduced Dirichlet boundary condition for the pressure at the inlet and outlet is proven to work very accurate, which is expressed in an error in the simulated pressure field of less than $0.002\,\%$ across the domain. The good agreement of the simulated velocity profiles (Fig. 5a) result from an accurate pressure field and valid formulations of the viscous and pressure forces (see Eqs. 5 and 6). The mean error of the local velocity at time $t_3^*$ is given with $\text {err}_u = 0.50\,\%$.

Fig. 5

a The analytic (lines) and computed velocity profiles (symbols) of the Poiseuille flow between the two plates for three time steps ($t^*_1 = 0.67$, $t^*_2 = 1.46$ and $t^*_3 = 4.17$. b The error of the computed pressure field over the length of the plates at $t^*_3$

Fig. 6

a The setup of the reduced model introducing all parameters is shown and b the combined momentum balance of the two liquids is given

4.1.2 Capillary Rise

In this simulation, the capillary rise between two flat parallel walls is examined. The height of the system is $H = 9.3\,\hbox {mm}$, and the gap between the plates is $d = 1.55\,\hbox {mm}$ wide. We compare our SPH simulation with a reduced model. The reduced model uses a global momentum balance of the two-phase system between the plates. The same approach but in a round capillary with only one fluid was used by Fries and Dreyer (2008). The model setup is shown in Fig. 6. The two main assumptions in the reduced model are a parabolic flow profile over the whole domain and a constant curvature at the interface. With the relation $h_n = H - h_w$ the two remaining independent variables are the height of the wetting phase $h_w$ and the contact angle $\theta$. Blake and Haynes (1969) proposed a linear correlation of the dynamic contact angle $\theta$ and the contact line velocity U in the vicinity of the static contact angle $\theta _0$:

$$\begin{aligned} U = \frac{\sigma _{wn}}{\xi }\left( \cos (\theta _0)-\cos (\theta )\right) . \end{aligned}$$

(24)

Here $\xi$ is called the intrinsic friction coefficient and has the unit $\hbox {Pa}\,\hbox {s}$. $\theta _0$ is the static contact angle and $\theta$ is the apparent/dynamic contact angle. This linear relation has also been used in our SPH simulation (Fig. 7a) with an intrinsic friction coefficient of $\xi = 0.041\,\text {Pa s}$ for the water–Fluorinert system. Blake (2006) found a similar friction coefficient for a comparable system ($16\,\%$ glycerol/water on PET). If the curvature of the interface is quasi-static during the capillary rise the assumption is valid, that the interface velocity $\mathrm{d}h_w/\mathrm{d}t$ and the contact line velocity U are the same. With this equation one can solve the reduced model with the information of the dynamic contact angle. To stay close to the respective experiment, we chose water ($\rho _n = 1000\,\hbox {kg}$, $\mu _n = 1.0\times 10^{-3}\,\hbox {Pa}\,\hbox {s}$) as light non-wetting phase and Fluorinert ($\rho _w = 1800\,\hbox {kg}$, $\mu _w = 4.7\times 10^{-3}\,\hbox {Pa}\,\hbox {s}$) as denser wetting phase. The static contact angle is the same as the one in the experiment $\theta _0 = 45^{\circ }$, and the surface tension is $\sigma _{wn} = 0.055\, \hbox {N}/\hbox {m}$, which corresponds to the surface tension between water and Fluorinert. The pressure boundary conditions where chosen to be $p_\mathrm{top} = 0\,\text {Pa}$ and $p_\mathrm{bottom} = \rho _n\,g\,H = 91.2\,\text {Pa}$. In Fig. 7 the good correlation of the evolution of the capillary rise between the reduced model and the SPH simulation is shown.

Fig. 7

a The simulated dynamic contact angles over the corresponding contact line velocities during the capillary rise and b the average level of Fluorinert during the capillary rise is plotted over time

4.2 Simulation of a Drainage Processes in a Porous Structure

The geometry of the porous structure in the simulation is the same as the one of the PDMS micro-model in the experiment (comp. Fig. 1). The length of the simulation domain including the inflow channel is $l=3.48\,\hbox {mm}$ with a resolution of the whole porous structure of $1200\times 362$ particles. So each particle represents an area of $2.9\times 2.9\,\upmu \hbox {m}$. We perform two different types of simulations. The first type is a dynamic drainage process with a fixed pressure difference between inlet and outlet. The other type is a quasistatic drainage process, where the displacement of the wetting phase is performed in multiple steps with relaxation times in between to measure system properties under static conditions.

4.2.1 Dimension Reduction

All simulations of the micro-model were performed in 2D in order to limit the computational load. By reducing the vertical dimension, we assumed a parabolic flow profile in vertical direction, to take the additional viscous friction into account and a constant curvature of the interface in vertical direction.

Viscous Force

The cross-sections of all capillaries in the micro-model have the shape of a rectangular duct. The viscous drag resulting from the no-slip condition at the bottom and top wall of the micro-model strongly influences the permeability $\kappa$ of the system. With the assumption of a parabolic flow profile in vertical direction, the flow of a rectangular duct is well approximated if the width b of the duct is larger than the height a. For this reason an additional viscous force is applied to take care of the dimension reduction in the simulation:

$$\begin{aligned} \frac{\mathrm{D}\mathbf {v}_{i,\mathrm{red}}}{\mathrm{D}t} = \frac{\mathrm{D}\mathbf {v}_{i}}{\mathrm{D}t} - \frac{12\,\mu _i\,\mathbf {v}_i}{\rho _i*a^2} \end{aligned}$$

(25)

The validation of this approach was carried out in a channel with a ratio $b/l = 3$. The height was defined as in the experiment $a = 100\,\upmu \hbox {m}$ and the pressure difference was chosen $\Delta p = -10\,\hbox {Pa}$. The width varied from $80\,\upmu \hbox {m} < b < 500\,\upmu \hbox {m}$, which corresponds to the minimum and maximum width of the channels in the micro-model. The average channel width of the micro-model was $b = 160\,\upmu \hbox {m}$. In Table 1 the influence of the dimension reduction is shown. The error in permeability of the mean pore size is $\approx$6 %, which is not negligible, but it decreases quickly in larger channels where the most dynamics take place. We included the comparison of the permeability between two flat plates (see last row in Table 1) to highlight the contribution of the dimension reduction in the other cases. The analytic solution for a flow in a rectangular duct is given for example by Tuller and Or (2001):

$$\begin{aligned} \begin{aligned} \bar{v}&= \frac{a^2\,B_s}{4\,\mu }\left( -\frac{\mathrm{d}p}{\mathrm{d}z}\right) \\ B_s&= \frac{1}{3} - \frac{a}{b} \frac{64}{\pi ^5} \sum _{n=1}^{\infty }\frac{\tanh \left( \frac{b}{a} \frac{\pi (2n-1)}{2} \right) }{(2n-1)^5}, \end{aligned} \end{aligned}$$

(26)

where we stopped the series at $n=1000$.

Table 1 Comparison of theoretic values of the permeability of rectangular ducts and simulation results using the dimension reduction

Surface Force

The capillary pressure of a rectangular duct is equal to

$$\begin{aligned} p^c = p^n-p^c=\frac{\sigma _{wn}}{r_c}, \end{aligned}$$

(27)

where the radius of the curvature $r_c = \frac{1}{\kappa }$ can be determined as (Joekar-Niasar et al. 2009):

$$\begin{aligned} \frac{8(\pi /4-\theta )r_c+[2(a+b)-8\sqrt{2}\cos (\pi /4+\theta )r_c]\cos (\theta )}{ab-4r_c^2[\sqrt{2}\cos (\pi /4+\theta )\cos (\theta )-(\pi /4-\theta )]} = \frac{1}{r_c}, \end{aligned}$$

(28)

where a and b is the height and the width of the channel and $\theta$ is the static contact angle. The surface tension between Fluorinert and water is $\sigma _{wn} = 0.055 \,\hbox {N}/\hbox {m}$. The static contact angle in the micro-model was identified as $\theta \approx 45^{\circ }$. In this case Eq. 28 simplifies to:

$$\begin{aligned} \kappa = \frac{1}{r_c} = \left( \frac{2}{a}+\frac{2}{b}\right) \cos (\theta ). \end{aligned}$$

(29)

Because the total curvature now adds up by the horizontal and vertical extension of the duct, the only necessary assumption for a reduction of dimension is a static contact angle of the interface on the top and bottom walls. Then we can get $p_c$ by computing only the horizontal curvature in the simulation and adding the constant contribution from the vertical curvature.

$$\begin{aligned} p_c = \sigma \kappa _h + 2\,\frac{0.055\,\hbox {N}/\hbox {m}\,\cos (45^{\circ })}{100\,\upmu \hbox {m}} = \sigma \kappa _h + 778\,\text {Pa} \end{aligned}$$

(30)

Fig. 8

a, b The phase evolution after $94\,\hbox {ms}$ with the corresponding dynamic pressure field during a drainage process with $1.86\,\text {kPa}$. c All wn interfaces are shown and d the experimental drainage process is shown after $1.45\,\hbox {s}$

4.2.2 Dynamic Drainage Process

For the dynamic drainage process, we performed several simulations with the pressure difference between inlet and outlet $\Delta p$ ranging from 1.678 to $5.778\,\text {kPa}$. Then we compared the simulated and experimental drainage process with $\Delta p = 1.86\,\text {kPa}$.

In Fig. 8 the simulated flow path after $94\,\hbox {ms}$ is shown together with the flow path in the experiment after $1.45\,\hbox {s}$. In Fig. 9 multiple simulated drainage processes with different applied pressure gradients are shown. We study the propagation of the interfacial area between wetting and non-wetting phase $A_{wn}$, the interfacial area between wetting and solid phase $A_{ws}$, the Helmholtz free energy F and the saturation of the wetting phase $S_w$. In an isothermal and isochoric system the change of the Helmholtz free energy, F is only affected by interface formation:

$$\begin{aligned} F(t) = A_{wn}(t)\,\sigma _{wn} - A_{ws,}(t)\,(\sigma _{ns}-\sigma _{ws}) + \mathrm{const}. \end{aligned}$$

(31)

Here $A_{wn}$ is approximately computed as product of the measured horizontal length of the wn interfaces in the 2D simulations

$$\begin{aligned} L_{wn,h} = \sum _i |\mathbf {n}_{wn,i}|\,V_0, \end{aligned}$$

(32)

where $V_0$ is the volume of one particle in the 2D simulation, and the constant vertical length of the curved interface

$$\begin{aligned} L_{wn,v}&= \frac{a\,(\pi -2\,\theta )}{2\,\cos (\theta )} = 111\,\upmu \hbox {m} \end{aligned}$$

(33)

$$\begin{aligned} A_{wn}&= L_{wn,h} \, L_{wn,v}. \end{aligned}$$

(34)

$A_{ws}$ is computed as sum of the vertical and horizontal ws-interfaces:

$$\begin{aligned} A_{ws}&= \sum _i \left( |\mathbf {n}_{ws,i}|\,V_0 \, a + 2\,\delta _{w,i} V_0 \right) . \end{aligned}$$

(35)

where $\delta _{w,i}$ is defined as:

$$\begin{aligned} \delta _{w,i} = {\left\{ \begin{array}{ll} 1 \quad \text {if particle i belongs to the wetting phase} \\ 0 \quad \text {if particle i belongs to the non-wetting phase} \end{array}\right. } \end{aligned}$$

(36)

During the dynamic simulation, we determined the capillary pressure from the output of the pressure field (see Fig. 8b). The capillary pressure was identified from the pressure difference at interfaces between the continuous wetting and non-wetting phase. Thus intrusions of the wetting phase had no effect on the measured capillary pressure.

Fig. 9

The different behavior of the dynamic drainage processes when applying pressure gradients between 1678 and $5778\,\text {Pa}$. We stopped the simulation in the case of $\Delta p = 1678\,\text {Pa}$ after the barrier was reached. In the other cases the end of the plotted line is equal to a breakthrough

Fig. 10

Dynamic drainage process: The figures in the left column show the flow path in the experiment over time. The images in the right column show the flow path in the simulation with $\Delta p = 1.86\,\text {kPa}$ at corresponding saturations

4.2.3 Quasistatic Drainage Process

The quasistatic drainage process in the simulation is setup as follows. The system is initially fully wetted with Fluorinert. The non-wetting phase is pushed with a certain velocity in the system. The boundary conditions are still Dirichlet boundary conditions for the pressure and homogeneous Neumann boundary conditions for the velocity as described in Sect. 3.4.2. We use a PI controller for the inlet pressure to obtain an average velocity in the inflow channel $\bar{v}_\mathrm{dyn}$. The setpoint of the average inflow velocity is set to $\bar{v}_\mathrm{dyn} = 1 \times 10^{-2}\,\hbox {m}/\hbox {s}$. The pressure at the outlet is fixed at $p=1.0\,\text {kPa}$.

The Fluorinert is displaced until the saturation $S_w$ has dropped by $5\,\%$ in the domain. After that the setpoint of the velocity is set to $\bar{v}_\mathrm{stat} = 0\,\hbox {m}/\hbox {s}$ for a certain time interval $\Delta t_\mathrm{stat} = 1\times 10^{-2}\,\hbox {s}$ to let the system relax. The capillary pressure was measured from the pressure difference between inlet and outlet in the relaxed system. Then the setpoint of the average inlet velocity is switched back to $\bar{v}_\mathrm{dyn}$.

5 Discussion of Results

We have performed both simulations and experiments in a geometric identical micro-model with the goal to get a detailed insight into the underlying physical processes. The CSF/CLF model in the SPH simulations was able to predict a flow path in the dynamic drainage process which is very close to the one found in the experiment (Fig. 10). The flow path in the simulated quasistatic drainage process is also similar to the one found in the experiment. Here one more pore body gets drained in the simulation when the saturation lowered to $S_w = 0.35$ (Fig. 11). The measured capillary pressures at specific saturations fit very well with the experiment for both the dynamic (Fig. 13a) and the quasistatic (Figs. 12, 14a) case. In both cases the peak in the capillary pressure at $S_w \approx 0.7$ is captured in simulation and experiment.

Fig. 11

Quasistatic drainage process: The figures in the left column show the system at rest under equilibrium conditions in the experiment. The images in the right column show static conditions in the simulation at corresponding saturations

Fig. 12

Comparison of the final states of the simulated drainage processes with $\Delta p = 1.86\hbox {\,kPa}$ and $\Delta p = 2.778\,\hbox {kPa}$

The SPH simulation can give very interesting insights into the dynamic drainage process (see Fig. 8). The variation of the local capillary pressure at the different interfaces can be monitored in detail. Also the interfacial area in 2D can be determined with an easy summation over the SPH particles. Nonetheless, the measured interfacial areas show a significant difference in both the dynamic and quasistatic drainage process. The comparison of the interfacial areas are shown in Fig. 13b for the dynamic drainage and in Fig. 14b for the quasistatic process. The main part of this difference results from the larger number of small intrusions of wetting fluid captured in the SPH simulations. Another part of the difference could also be explained by the different evaluation processes. In the simulation the interfacial area is computed by a simple product of the measured interface in 2D and a constant contribution of the height of the micro-model, which will overestimate the interfacial area by some percent.

In the simulations, we performed several drainage processes with increasing pressure difference (Fig. 9). Apart from a faster decrease in the saturation of the wetting phase, also the total change in saturation at breakthrough varies. In Fig. 12 a comparison between the final states before breakthrough of the simulated drainage process with an applied pressure difference of $\Delta p = 1.86\,\text {kPa}$ and $\Delta p = 2.778\,\text {kPa}$ is shown. With a pressure difference of only $\Delta p = 1.86\,\text {kPa}$ one can see a clear effect of capillary fingering, both in the simulation and the experiment (compare final states in Fig. 10). This effect vanishes at pressure gradients of $\Delta p = 2.778\,\text {kPa}$ and higher. This observation corresponds to the work of Joekar-Niasar (2012), who showed in a non-equilibrium pore network simulation that pressure–saturation curves depend on the capillary number.

From Fig. 9 one can also find a remarkable difference in the interfacial areas between the cases of $\Delta p = 2.278$ and $5.778\,\text {kPa}$. Even the amount of displaced wetting phase is very similar, and $A_{wn}$ gets approximately twice as large in the case of $\Delta p = 5.778\,\text {kPa}$, while $A_{ws}$ stays smaller. The reason for this effect is remaining films of wetting phase in pores which have been drained very fast.

Fig. 13

a A capillary pressure–saturation plot for the dynamic drainage process is shown for both simulation and experiment. b The evolution of the measured interfacial area between wetting and non-wetting phase is plotted

Fig. 14

a The capillary pressure–saturation plot for the quasistatic drainage process is shown for both simulation and experiment. b The evolution of the measured interfacial area between wetting and non-wetting phase is plotted

Fig. 15

a The propagation of the saturation during the dynamic drainage process in the experiment and simulation. b The propagation of the simulated drainage process is shown in detail

A major disagreement between simulation and experiment is found in the dynamics of the displacement process (Fig. 15). We assume two different effects, we have not captured in the simulation, causing the higher volume flux in the simulation. The first effect is an additional flow resistance in the experimental micro-model setup, which could originate for instance from a flow resistance at the connections of the micro-model with the reservoirs. To capture this resistance, we compared the measured and simulated volume flux for a single-phase flow of dyed water. With a fixed pressure gradient of $\Delta p = 1.86\,\text {kPa}$ the simulated volume flux was $0.88\, \hbox {ml}/\hbox {min}$ where the measured volume flux was only $0.32\, \hbox {ml}/\hbox {min}$. So the permeability of the whole micro-model setup in the experiment is about 2.75 times lower than in the simulation. With this knowledge, we can take the simulated saturation profile and add an additional resistance term for single-phase flow (RSP) which is linear to the current reciprocal volume flux (cf. Fig. 15b):

$$\begin{aligned} \frac{1}{\dot{V}_\mathrm{RSP}} = \frac{1}{\dot{V}_\mathrm{SIM}} + a \frac{1}{\dot{V}_\mathrm{SIM}}, \end{aligned}$$

(37)

with $a = 1.75$. While the dynamics of the saturation change for the saturation region of $S_w < 0.6$ now become similar with the additional resistance, this is not true in the time before (‘RSP’ in Fig. 15a).

In the experiment the bottle effect at a saturation of $S_w \approx 0.7$ is holding the process at almost rest for one second, where the effect is by far not that strong in the simulation. The chosen pressure difference for the comparison with $\Delta p = 1.86\,\text {kPa}$ lays just slightly above the theoretic entry pressure at $S_w \approx 0.7$. Therefore we expect another effect called “stick–slip” (Blake and Shikhmurzaev 2002), which is associated with slow contact line movement, being responsible for the remaining offset. The “stick–slip” region is characterized by a strong dependence of the dynamic contact angle on the capillary number (Ca). While the CLF model properly captures the dynamic contact angle at high Ca-numbers, it can’t capture the “stick–slip” behavior related to contact angle hysteresis on rough surfaces. To get an idea how this “stick–slip” behavior could change the drainage process a second resistance is added. This resistance at the contact line (RCL) increases with lower volume fluxes. The combined resistance of single-phase flow and the contact line could be expressed by a power law of the form:

$$\begin{aligned} \frac{1}{\dot{V}_\mathrm{RSP \,{+}\,RCL}} = \frac{1}{\dot{V}_\mathrm{RSP}} + b \left( \frac{1}{\dot{V}_\mathrm{RSP}}\right) ^{c}. \end{aligned}$$

(38)

With $b= 7.8 \times 10^{-3}$ and $c=2.5$ the modified drainage process gets very similar to the experimental one (‘RSP $+$ RCL’ in Fig. 15a). Therefore, a more detailed analysis of the wetting dynamics of water and Fluorinert on PDMS for low Ca-numbers is necessary for a better understanding of the “stick–slip” behavior of the system of interest. This study illustrates that a detailed knowledge of the dynamic contact angle is also important when simulating flow with low Ca numbers, where static contact angles can be assumed constant most of the time. With a model to describe the dynamics of the contact line separated from the interface, like the CLF model, it is possible to include information on roughness in an extended model description in the future.

6 Conclusion

We performed simulations and experiments of a drainage process in a small micro-model. The implemented SPH model is able to describe the effects of surface tension in a multi-phase system. We found good agreement between simulations and experiments for the flow path and the evolution of the capillary pressure. Combined with the easy tracking of interfacial areas, SPH simulations are an appropriate tool to capture $p_c$–$S_w$–$a_{wn}$ relationships for larger geometries in the future. The uniqueness of those relationships could be studied during drainage and imbibition, which was proposed by Hassanizadeh and Gray (1993). Apart from those good correlations a major difference in the dynamic evolution of the drainage process was found. We assume surface roughness and sticking effects play a major role when displacing a wetting fluid from a porous structure. These effects need further investigations.