Assessing the Mobilization of Trapped Mass of Emulsions Flowing in an Idealized Pore Using the Lattice Boltzmann Method

Abstract

In this work, a Lattice Boltzmann model for multi-component fluids augmented with repulsive near-contact interactions is employed to simulate the dynamics of flowing emulsions within idealized pores. The model is firstly tested on experimental data of droplets’ self-assembly in diverging-converging micro-channels from literature and then used to investigate the trapping/mobilization of the dispersed phase of an emulsion in an idealized series of pores, as influenced by both the Capillary number and the solid walls wettability. Both parameters may vary as the result of an injection of surfactants, a procedure commonly adopted in soil remediation and Enhanced Oil Recovery applications. The analysis shows that the proposed model is able to reproduce correctly the experimental data and gives interesting insights on the trapping/mobilization phenomenon resulting from a modification of the flow conditions caused by the injection of surfactants.

Article Highlights

• A LBM multiphase model is employed to analyze a two-phase flow in a simplified pore geometry.

• Capillary number and solid wettability are varied to simulate the injection of surfactants.

• The mobilization of the previously trapped mass is quantified and related to the modified flow conditions.

1 Introduction

A two-phase emulsion is a fluid system consisting of two interacting immiscible fluids: one, referred to as the dispersed phase, is composed by droplets of various sizes carried by the second, called continuous or carrier phase. Such systems play a pivotal role in many areas of environmental fluid mechanics, from petroleum and biochemical engineering to contaminant hydrology (Perrazzo et al. 2018). In particular, it has been proven that the degree of emulsification may be increased and/or stabilized by means of surfactants (chemicals that affect the capillary characteristics of the solid–liquid interface), to improve the efficiency of processes like the Enhance Oil Recovery (EOR) (Hirasaki et al. 2011; Muggeridge et al. 2014), as well as aquifers remediation techniques (Harwell et al. 1999; Javanbakht et al. 2017).

These procedures are aimed at extracting the fraction of dispersed phase which adheres to the walls of the porous matrix. The trapping phenomenon depends on both the geometry of the pores and the flow conditions (Perrazzo et al. 2018; Shams et al. 2021). The latter comprises the contact angle with the solid walls and surface tension, both being affected by the presence of surfactants (Wasan and Nikolov 2003; Kralchevsky et al. 2015). Indeed, the composition of the emulsions-surfactant solution affects the surface tension, which is connected to the contact angle according to Young’s equation (Reed and Healy 1984; de Gennes 1985). However, the thorough comprehension of the emulsions dynamics in porous matrices, and its quantitative response to the presence of surfactants, in terms of the relevant physical parameters, is still an open issue in the design of soil remediation campaigns (Yang et al. 2022).

Flooding tests have been widely employed to analyze the effect of surfactants on oil displacement (Dwarakanath et al. 1999; Lu and Pope 2016); however, such studies mainly focused on macroscopic parameters (such as pressure, displacement efficiency, etc.), while effects of the surfactant phase on microscopic displacement dynamics were missing (Yang et al. 2022). The adoption of micro-CT technology within the flooding tests allowed to evaluate important parameters for the assessment of the Capillary number (Georgiadis et al. 2013; Armstrong et al. 2014; Zhu and Wang 2017), the latter expressing the ratio between viscous and surface tension forces (Melrose and Brandner 1974; Hilfer and Øren 1996). However, such technology requires expensive and precise equipment (Wang et al. 2019) and is not-fully reliable (Garfi et al. 2020; Huang et al. 2021). Computational modelling of flowing emulsions in porous media can make up for the drawbacks of experimental approaches. Several models have been proposed featuring different levels of complexity.

The least complex numerical modelling is the Pore Network Modelling (PNM) (Blunt and King 1991), where both the porous matrix and the emulsion hydrodynamics are strongly simplified: in particular, the former is assumed to be represented by a network of channels and cavities of regular shapes (Dias and Payatakes 1986; Mogensen and Stenbt 1998; Dong and Blunt 2009). These models allowed to analyze fundamental aspects of emulsions dynamics (Valvatne and Blunt 2004; Niasar et al. 2009; Wang et al. 2015); nonetheless, it has been found out that the over-simplifying assumptions may lead to several inaccuracies (Niasar and Hassanizadeh 2012; Wei et al. 2020).

More sophisticated models, based on Navier–Stokes equations in complex domains (Wörner 2012), were able to overcome such problems, allowing for a more detailed description of the phases and solid interactions. However, they had to be coupled with interface tracking algorithms (Jettestuen et al. 2013; Ferrari and Lunati 2013; Friis et al. 2019), often resulting in a remarkable computational burden (Wang et al. 2019).

Recently, multi-component versions of the Lattice Boltzmann Method (LBM) (Succi 2001; Krüger et al. 2017; Wei et al. 2018) has been extended to deal with various phenomena connected to multiphase flows (Radhakrishnan et al. 2022), from the effects of the Capillary number on fluid displacement in natural soil (Tsuji et al. 2016) and the estimation of relative permeability (Jiang and Tsuji 2017) to the effect of different wettability conditions (Wei et al. 2019, 2020; Mora et al. 2021b, a) and surfactants dynamics (Liu et al. 2018; Zhu et al. 2019; Zhang et al. 2021), paving the road for the investigation of emulsion flows at the pore scale. Furthermore, a numerical method like LBM, thanks to its mesoscopic approach, stands as a computationally efficient alternative to standard Navier–Stokes solvers due to its coding simplicity and amenability to parallellization (Harting et al. 2005; Bernaschi et al. 2011; Succi et al. 2019).

Hence, in this study, the multi-component LBM model proposed in Montessori et al. (2019), incorporating near-contact interactions, is adopted to investigate emulsion dynamics in an idealized porous matrix. This work aims at assessing the role of the Capillary number and solid walls wettability on the mobilization process of trapped dispersed phase within simplified pores; these parameters are varied in order to mimic the injection of surfactants in the continuous phase of the emulsion. The analysis is also aimed at understanding if this class of models may be a valuable tool for the evaluation of operational scenarios in soil remediation and oil extraction.

The paper is structured as follows: firstly, a brief description of the adopted LBM model is presented; secondly, the model is further validated against experimental data from literature; thirdly, the investigation of the trapping process is carried out. Finally, some conclusions are drawn.

2 Method

In this work, the multi-component LBM model, based on the color-gradient approach (Leclaire et al. 2017) and developed by Montessori et al. (2019) (to which the reader is addressed for a detailed description) is employed; exhaustive reviews on the LBM can be found in Succi (2001), Sukop and Thorne (2006) and Krüger et al. (2017). Some additional information about the adopted model can be also found in Appendix A.

An important feature of the adopted model is the inclusion of a short-range repulsive forcing term, able to mimic Near-Contact Interactions (NCIs); such NCIs may arise as an ultimate consequence of different conditions (such as ions absorption, electric charges or the presence of chemicals on interface surfaces) and are needed to withstand the effect of the surface tension by frustrating the coalescence between interacting interfaces.

The two emulsion’s fluid components are described by two distinct sets of Probability Distribution Functions (PDFs) $f_a^k(\textbf{x},t)$, $k=1,2$, which represent the probability of finding a fluid particle of the $k^{th}$ component at position $\textbf{x}$, with velocity $\textbf{e}_a$ at time t; a is an index spanning over the lattice discrete directions. In this case a D3Q27 lattice is adopted (Krüger et al. 2017), thus $a=0,...,26$.

The PDFs are governed by the equations:

$$\begin{aligned} f_a^k(\textbf{x}+\mathbf {e_a}\Delta t, t+\Delta t) = f_a^k(\textbf{x},t)+\Omega _a^k[f_a^k(\textbf{x},t)] \end{aligned}$$

(1)

The density of the $k^{th}$ component and the total momentum of the emulsion is given by the zero$^{th}$ and first order momentum of the PDFs as:

$$\begin{aligned} \rho ^k(\textbf{x},t)=\sum _af_a^k(\textbf{x},t) \qquad \rho (\textbf{x},t)\textbf{u}(\textbf{x},t)=\sum _k\sum _a f_a^k(\textbf{x},t)\textbf{e}_a \end{aligned}$$

(2)

The term $\Omega _a^k[f_a^k(\textbf{x},t)]$ is the collision operator and is defined as:

$$\begin{aligned} \Omega _a^k=(\Omega _a^k)^{(3)}\Big [(\Omega _a^k)^{(1)}+(\Omega _a^k)^{(2)}\Big ] \end{aligned}$$

(3)

$(\Omega _a^k)^{(1)}$ is the standard relaxation operator and is connected to the kinematic viscosity $\nu$ of the $k^{th}$ component (Succi 2001). The term $(\Omega _a^k)^{(2)}$ is the perturbation operator and accounts for the interfacial tension $\sigma$ (Gunstensen et al. 1991). Lastly, $(\Omega _a^k)^{(3)}$ is the so-called recolouring operator, regulating the mutual diffusion of the components (Gunstensen et al. 1991).

It is possible to obtain the Navier Stokes equation from Eq. (1) by performing the Chapman-Enskog expansion (Chen and Doolen 1998).

The stress jump across a fluid interface can be calculated as:

$$\begin{aligned} \textbf{T}^1\cdot \textbf{n}-\textbf{T}^2\cdot \textbf{n}=-\nabla \cdot (\sigma \textbf{I}-\sigma (\textbf{n}\otimes \textbf{n}))-\pi \textbf{n} \end{aligned}$$

(4)

where $\textbf{T}^k=-p^k\textbf{I}+\rho ^k\nu (\nabla u+\nabla u^T)$ is the stress tensor of the $k^{th}$ fluid component and $\textbf{n}=-\nabla \psi /|\nabla \psi |$ is the unit vector normal to the interfacial surface. The scalar $\psi$ is the so-called phase field, defined as $\psi =(\rho ^1-\rho ^2)/(\rho ^1+\rho ^2)$. The last term $\pi [d(x)]$ represents the repulsive near-contact force acting on the interface; d(x) being the distance between two different interfaces along the surface normal direction $\textbf{n}$.

Assuming $\textbf{T}\simeq -p\textbf{I}$ and projecting Eq. (4) along the normal, the augmented Young-Laplace equation is obtained (Chan et al. 2011):

$$\begin{aligned} p^1-p^2=\sigma (\nabla \cdot \textbf{n})-\pi \end{aligned}$$

(5)

The term $\pi$ can be included in the LBM framework considering a near-contact repulsive forcing term, which acts only on the interface, defined as:

$$\begin{aligned} F_{rep}=\nabla \pi :=-A_h[d(\textbf{x})]\textbf{n}\delta _I \end{aligned}$$

(6)

being $\delta _I=\frac{1}{2}|\nabla \psi |$ the function that confines the forcing term effects only on the interface and $-A_h[d(\textbf{x})]$ the magnitude of such NCI (a force per unit of volume); $A_h$ is a function of the distance $d(\textbf{x})$ and is a constant if $d<d_{min}$, while it vanishes for $d>d_{max}$. The force term in Eq. (6) can be included in Eq. (1) as usual (Krüger et al. 2017).

For the NCIs approach an important dimensionless parameter is the ratio of near-contact versus capillarity forces $N_c$, defined as: $N_c=A_h\Delta x/\sigma$, with $\Delta x$ the grid spacing and $\sigma$ the surface tension (see Montessori et al. (2019, 2020, 2021) for further details). In this work, we set $N_c \sim O(0.1\div 1)$, which is comparable with the expected physical values (Israelachvili and Pashley 1982; Donaldson et al. 2011).

Finally, as shown in Latva-Kokko and Rothman (2005b), in the color-gradient LBM approach it is possible to set the cosine of the contact angle $\vartheta$ in order to take into account the solid wall wettability. Indeed, the algorithm simulates the presence of a fictious liquid phase by setting in the solid nodes the density $\rho ^k$ of the $k^{th}$ component to an artificial value; the variation of $\vartheta$ between 0$^\circ$ and 180$^\circ$ is achieved by imposing $0\le \rho ^k\le 1$. A remarkable advantage of this approach is that $\vartheta$ can be varied independently from other quantities.

3 Validation: Droplets Arrangement Inside a Pore

The validity of the proposed model, which has been already adopted to analyze important features of multiphase flows (Tiribocchi et al. 2019; Montessori et al. 2020, 2021), is here assessed against the experimental data from Jose and Cubaud (2012); the aim is to numerically replicate the droplets arrangements observed in the experiments.

The experimental set-up of Jose and Cubaud (2012) is composed by three consecutive elements:

• 1. a square inlet micro-channel;

• 2. a prismatic microfluidic chamber;

• 3. an outlet square channel;

Fig. 1

Computational domain in the validation case: in blue a sketch of the sections of the internal periodic conditions

The experiments consist in the analysis of the dynamics of droplets inside the prismatic chamber in different flow conditions. Droplets feature an initial diameter $D_0$ which is comparable to the inlet/outlet channel width h and are generated by means of a T-junction located at the beginning of the inlet channel. The viscosity of the dispersed phase, i.e. the droplets, is $\mu _o$, while the viscosity of the carrier phase is $\mu _w$. Different pairs of fluids with different viscosity ratios $\chi =\mu _o/\mu _w$ ($1.9\le \chi \le 44$) are used. An important parameter regulating the overall dynamics is the Capillary number Ca, defined in terms of the viscosity of the carrier fluid $\mu _w$ as in Jose and Cubaud (2012):

$$\begin{aligned} Ca=\frac{U_c\mu _w}{\sigma } \end{aligned}$$

(7)

where $U_c$ is the mean velocity inside the inlet channel. Ca represents the ratio of the orders of magnitude of the viscous and surface tension force. In the case of the experiment of Jose and Cubaud (2012) the walls of the device are not wetted by the dispersed phase.

Furthermore, the Reynolds number of the carrier fluid inside the pore chamber $Re_w=U_pD_0\rho _w/\mu _w$ (with $\rho _w$ the carrier density and $U_p$ the mean velocity in the middle of the pore) is small ($Re_w<<1$).

3.1 Simulation Set-up

The numerical domain is shown in Fig. 1: the channel width h is set to 20 lattice units (lu in the following), being 1 lu equal to 12.5 $\mu$m; the thickness of the domain, which is not represented in the figure, is set to 20 lu.

The inlet channel is equipped with an internal periodic boundary condition which is employed to generate droplets trains featuring the desired spacing; such condition is obtained by copying the PDFs of the two phase from the internal section to the upstream boundary (see Fig. 1). The outflow at the outlet channel (right) is implemented through an open boundary condition (Krüger et al. 2017); no-slip conditions are applied on solid nodes and other boundaries (Succi 2001).

Table 1 Simulation parameters in the LBM model validation

In order to replicate the experimental condition a Capillary scaling is assumed ($Ca_{num}=Ca_{exp}$) as well as a matching of the viscosity ratio. The numerical viscosity ratio $\chi _{num}$ is set for all the simulation equal to 1.9 (with $\mu _{w,num}$ assumed to be equal to 0.167). The numerical surface tension $\sigma _{num}$ is set to 0.05 which is compatible with the dynamics under investigation. The flow is obtained by applying on both phases a mass force, which is tuned to obtain a $U_{c,num}$, evaluated a posteriori, that satisfies the Capillary scaling. The other parameter controlling the droplets arrangements inside the chamber is the ratio $L_0/D_0$, where $L_0$ is the initial distance between two droplets and $D_0$ is their diameter inside the inlet channel. The inflow numerical droplet diameter is set to $D_{0,num}\simeq 0.9h$, consistently with the experimental values; $L_0$ can then be imposed by adjusting the location of the internal periodic boundary condition (as shown in Fig. 1). The values of such parameters explored in the numerical simulations, along with the respective Reynolds numbers, are listed in Table 1.

Fig. 2

Phase diagram of the benchmark simulations. Circles: experimental droplet arrangements in Jose and Cubaud (2012). Stars: LBM model simulations. The color of both experimental and numerical marker depicts the attained droplet arrangement. The filled gray area is the Capillary number range considered in the analysis of mobilization in Sect. 4

In Fig. 2, the observed and simulated droplets arrangements are reported in the $Ca-L_0/D_0$ plane. Each type of arrangement is depicted by a different marker color; LBM simulations and experiments are represented by stars and circles respectively; as shown, simulations have been carried out spanning the range of Ca considered in the analysis discussed in the next section (roughly $0.02<Ca<0.2$) and represented by the filled gray area of the graph.

3.2 Benchmark Results

The phase-diagram depicted in Fig. 2 shows that the LBM model is able to reproduce the different experimental arrangements; in Fig. 3 a visual comparison between simulations and experiments is shown. It is worth specifying that spurious small droplets, which appear in the experiments and are a consequence of T-junction generation, are not reproduced by the LBM model since the generation process is not replicated.

Fig. 3

Comparison between physical results from Jose and Cubaud (2012) (panel a) and LBM simulations (panel b). Panel b: $\psi$ field is shown; the red boundary of the droplets shows the local thickness of the interface

The overall agreement is fairly good: a close match with the experimental results can be noticed for cases with one and two layers droplets. A more quantitative assessment is shown in Fig 4, where the variation of the distance between two droplets $L(\tilde{x},t)$ inside the computational domain for the one-ordered layer is depicted. For a given pair of droplets the distance is computed as $L(\tilde{x},t)=x_R(t)-x_F(t)$ with $x_R$ the rear position of the first droplet and $x_F$ the front position of the second; the spatial coordinate $\tilde{x}$ of each distance is then evaluated as their mean position ($\tilde{x}=(x_R(t)+x_F(t))/2)$. The LBM model is able to reproduce such spatial variation quite well.

Fig. 4

Spatial distribution of the distance L between two droplets for one-layer arrangement: data from Jose and Cubaud (2012) (red circles) are compared with LBM values (blue solid line)

The case of one ordered layer exposes a well-known, general weakness of the outflow boundary conditions for multiphase models (Lou et al. 2013): a spurious higher outflow rate for the dispersed phase which, in turn, induces the slight asymmetry of the distribution of distances, as shown in Fig. 4. The problem is still an open issue. For the case considered in this work no major influence other than the one observed above has been detected, and we chose the type of boundary condition (namely the anti-streaming) that minimizes upstream reflections and provides a parameter-free implementation.

As per the remaining arrangements in Fig. 3 (VAL5 and VAL3), it is worth highlighting that, as reported in Jose and Cubaud (2012) and verified numerically, the cases resulting in multiple droplets layers (more than three) show abrupt transitions between states induced by even small deviations from the nominal values of Ca and $L_0/D_0$. Given such chaotic transition from three to six layers, the quantitative assessment is carried out by measuring lumped quantities such as the spatial distribution of the droplets arrangement amplitude A (normalized over the diameter inside the pore D) as a function of the chamber length x/L for different ratio $L_0/D_0$.

The stream of droplets inside the pore at the steady state is contained within an envelope of amplitude A that increases with the increase of the concentration of the dispersed phase (Jose and Cubaud 2012), i.e. with the decrease of $L_0/D_0$, and whose shape is well approximated by the parabola:

$$\begin{aligned} \frac{A}{D}=4\frac{x}{L}(1-\frac{x}{L})\frac{A_{max}}{D} \end{aligned}$$

(8)

With $A_{max}$ the maximum amplitude of the envelope, L the length of the pore chamber ($L=19h$).

Figure 5a shows the dimensionless amplitude A/D achieved by the droplet envelope as a function of the dimensionless distance x/L from the pore inlet. The fact that the amplitude A of the envelope of droplets increases with the decrease of $L_0/D_0$ is illustrated in Fig. 5b. The comparison between the experimental values and the numerical results, shown in both panels of Fig. 5, is fairly good.

Fig. 5

Panel a: spatial distribution of droplets arrangement amplitude A/D: experimental values from Jose and Cubaud (2012) (lines) are compared with LBM results (lines with star markers). Maximum amplitude attained in the chamber $A_{max}/D$ as a function of $L_0/D_0$ (panel b): comparison between experimental values from Jose and Cubaud (2012) (circles) and LBM results (stars)

The above analysis confirms the LBM model ability to correctly reproduce complex droplet arrangements in a confined geometry.

4 Mobilization of Slowly Moving Mass of Dispersed Fluid in a Series of Idealized Pores

In this section the LBM model is employed to study the dynamics of an oil-in-water emulsion flowing through an idealized series of pores. The aim of this analysis is to quantify how the variations of the capillary number Ca and the wettability angle of the solid walls $\theta$, promote the mobilization of the trapped/slow moving mass fraction of the dispersed phase. Indeed a variation of Ca and $\theta$ may be regarded as the primary effect of an injection of surfactants.

4.1 Simulation Set-up

The numerical domain is reproduced in Fig. 6: a narrow channel (the pore throat) of width h ($h=20$ lu) and length 2h links two half pores having the same dimensions and shape of the previous section (the domain thickness is set to 20 lu). In the following, the visualization of the domain is centered at the pore throat, in order to focus on the location where the dispersed mass develops most of the relevant dynamics. Periodic boundary conditions are applied at the inlet (left) and outlet (right) in order to reproduce a series of idealized pores; no-slip conditions are applied on solid nodes and other boundaries.

Fig. 6

Computational domain and initial condition used in the LBM model: the initial drop (the dispersed phase) is coloured in yellow, the carrier phase in blue and the solid matrix in cyan

As per the initial conditions, shown in Fig. 6, the two immiscible fluids are confined in the two halves of the chamber: the yellow fluid (i.e. the dispersed phase with viscosity $\mu _o$) in the upstream half and the blue fluid (i.e. the carrier phase with viscosity $\mu _w$) in the downstream half. As in the validation case, the viscosity ratio $\chi$ is set to 1.9, the motion is driven by a constant and uniform mass force and the Reynolds number of the carrier inside the pore is small ($Re_w<<1$).

The dynamics of the emulsion is analyzed in terms of Ca and the wettability of the pore’s walls, quantified by the contact angle $\vartheta$. Several combinations of these two parameters are tested and listed in Table 2.

The numerical simulations are performed so that a steady state is achieved and stably maintained. Steady state is assumed to be reached once the adhered mass at the solid walls does not vary further.

Table 2 Values of $\vartheta$ and Ca adopted in the LBM simulations

4.2 Discussion

In order to assess the degree of mobilization attained at each combination of parameters, the following quantities related to the dispersed phase are defined: (i) the adhered mass $M_{ad}$ and the number of free flowing droplets $N_{fd}$, (ii) the trapped mass. The adhered mass is the one in contact with the solid walls, and, together with the free flowing droplets, forms the total dispersed mass $M_{tot}$. The trapped mass is here defined to be the fraction of the dispersed mass flowing with a velocity lower than a certain threshold.

4.2.1 Adhered Mass and Free Droplets

The percentage of adhered mass $M_{ad}/M_{tot}$ in the LBM simulations is displayed in Fig. 7; since in simulations 01, 04 and 06 no adhesion has been imposed ($\vartheta =180^\circ$), the percentage is equal to zero. As expected, cases characterized by low adhesion conditions (high values of Ca and $\vartheta$) correspondingly attain small values of adhered mass at the steady state. Nonetheless in case 02, with $Ca=0.0423, \vartheta =0^\circ$, the steady adhered mass is smaller than that of case 05, with $Ca=0.23,\vartheta =0^\circ$. This fact can be explained with reference to Fig. 8, which shows the number of free droplets $N_{fd}$ generated by the break-up of the initial mass (panel a) and their median radius normalized with the channel width $R^*_{med}$ (panel b) (simulations 01, 04 and 06 are not included since there is no adhesion). As shown, while the simulations 02 and 05 share a similar amount of droplets at steady state, in case 02 the droplets radius is bigger than in case 05 leading to a higher amount of free mass. Both an increase in velocity and a decrease in surface tension promote the break up of larger droplets into smaller ones. Therefore, the higher the Ca, the more free droplets are expected to be generated by the break-up of the initial mass, however, as far as the results of this study are concerned, this correlation seems to be loose.

Fig. 7

Percentage of adhered dispersed mass ($M_{ad}/M_{tot}$) at the stationary state in LBM simulations

Fig. 8

a Number of free droplets $N_{fd}$ generated from the initial drop at the stationary state; b normalized median droplet radius $R^*_{med}$ at the steady state

It is worth noticing that in some simulations, the break-up of the free dispersed mass into smaller droplets is more evident than in other cases (e.g. in cases 02 and 05). This fact may seem counter-intuitive, however it appears to depend on the interaction between the adhered film created at the pore throat and the free mass: the build up of a slow moving film of dispersed phase at the wall induces a reduction of the effective throat section which, in turn, promotes the premature break-up in the form of dripping in presence of NCIs. Further investigations, object of future studies, are nonetheless needed in order to fully understand this kind of event.

4.2.2 Trapped Mass

As previously mentioned, the adhered mass should not be regarded as the only proxy for trapped mass. Indeed, a much richer picture of the overall emulsion dynamics can be obtained by investigating the dispersed mass distribution in velocity at a gauging section x, as shown in Fig. 9. Such quantity provides a comprehensive description of how fast the dispersed phase is moving through the porous system. Each bar represents the average dispersed mass percentage $M^*_x(u^*)$ moving with a longitudinal non-dimensional velocity $u^*=u/u_{max}$ where $u_{max}$ is the maximum longitudinal velocity reached in the section. It is worth clarifying that $M^*_x(u^*)$ is evaluated by averaging the dispersed mass percentage of different simulation snapshots once the steady state is reached.

As expected, high adhesion conditions, characterized by low Ca and low $\vartheta$ (simulation 07), entail a large amount of almost still mass. On the other hand, no adhesion conditions, characterized by high Ca and high $\vartheta$ (simulation 04), call for high velocities of the dispersed phase. Overall, a shift of the peak of the histogram towards high velocities, that is a mobilization of the trapped mass, can result from either a rise of Ca or an increase of $\vartheta$, or from a combined effect.

Fig. 9

Dispersed mass distribution in velocity; $M^*_x$ is the mass percentage, $u^*$ is the dimensionless longitudinal velocity. Red bars represent the trapped dispersed mass (assumed to be the one moving with $u^*\le 0.05$). Panels are positioned according to increasing Ca (left to right) and contact angle $\vartheta$ of the dispersed phase with the solid walls (bottom to top). Each inset shows the distribution of the two phases (see the caption of Fig. 6 for the interpretation of colours) in the pore geometry at an instant representative of the reached steady state; the red vertical cross section is the location of the mass fraction measurement

Furthermore, it is possible to notice that such shifts of the peaks, either caused by a rise in Ca or by a increase in contact angle, is always associated with a decrease in peaks height: a more oleophobic solid matrix promotes faster sliding of adhered dispersed mass. Additionally, this mass distribution in velocity may be informative of the emulsion dynamics: for example, the presence of a bimodal peak (as in simulations 03 and 08) reveals an asymmetry of film adhesion across the pore; it has been demonstrated that such instabilities in multiphase flows happen even at low Re (Stolovicki et al. 2018; Datta et al. 2022) and in this case are likely to depend on both the neutral wetting, adopted in these simulations, and the break-up of the initial drop. Moreover, the existence of free flowing droplets, as the ones generated in simulations 02 and 05, is revealed by smaller secondary peaks to the right of the main one.

A comparison between the profiles of the longitudinal velocity $u^*$ at the gauging section x for case 07 and 04 is shown in Fig. 10; the difference between the two dynamics is quite remarkable: while case 07 is characterized by thick layers of slow moving dispersed phase close to the wall, in case 04 such layer is not present. This description explains the high percentage of trapped dispersed mass in simulation 07, represented by the red bars in Fig. 9 and here assumed to be the one moving with $u^*\le 0.05$ (i.e. one order of magnitude lower than the average velocity in the section). Such behavior stands in clear contrast with the outcome of simulation 04, where this percentage is negligible. In general, it is possible to state that low $\vartheta$ and Ca cause the formation of a layer of slow moving mass whose thickness decreases as the two parameters grow larger.

Regarding the threshold $u^*\le 0.05$, it is worth clarifying that it was chosen in order to be representative of the high residence times connected to the trapped dispersed mass. Besides, it has been evaluated that, at least for this case study, the threshold sensitivity is quite low, hence the usage of lower thresholds would have given the same results.

Fig. 10

Profile of the longitudinal velocity $u^*$ for case 07 (panel a), 04 (panel b) at the gauging section x

It is possible to summarize the effect of the two parameters on the trapped dispersed mass by constructing plots like the one shown in Fig. 11. The maximum value of trapped dispersed mass percentage (over $70\%$) is yielded by simulation 07, featuring the lowest $\vartheta$ and Ca. The influence of the Capillary number is remarkable: doubling the value of Ca roughly halves the trapped mass, as shown by comparing case 07 with case 02 and case 09 with case 03; moreover, a value ten times higher causes the mobilization of the whole dispersed mass (as in case of simulation 05). Furthermore, even with low values of Ca, full oleophobic walls, as in simulation 06, highly reduce the trapping phenomenon, with a percentage of moving mass of over $85\%$. Hence, despite the coarse spanning of parametric space carried out in this study, it is clear how the impact of the contact angle $\vartheta$ on the trapping phenomenon is remarkable. However, simulations results confirm that even a small variation of Ca promotes a marked mobilization of the slowly moving mass.

The results, reported so far and condensed in Fig. 11, show the capability of the proposed LBM model to quantify the trapping phenomenon. The generated synoptic maps, such as the one shown in Fig. 11, may serve as guidelines for the implementation of operational protocols for oil extraction or soil remediation. Further investigations should focus on the assessment of the influence of complex geometries/heterogeneity on the trapping/mobilization dynamics.

Fig. 11

3D bar plot of the percentage of the trapped dispersed mass inside the pore in the parametric space spanned by the simulations

5 Conclusion

Emulsions flowing in porous matrices play an important role in many engineering fields, including oil extraction and remediation of contaminated aquifers. Nowadays, operational protocols are mostly formulated based on results from previous experience and flooding tests, both approaches exposing several flaws, costs and scaling issues being the most prominent ones. Resorting to a numerical approach seems therefore beneficial.

In this study, an existing multiphase numerical model based on the multi-component Lattice Boltzmann Method (LBM) color-gradient model, capable of modelling the short-range repulsive forces that aid to hinder the coalescence between interacting interfaces, is applied to study the trapping phenomenon in an idealized porous matrix.

Firstly, the model is validated against experimental data concerning droplet arrangements inside a pore. Then, it is used to quantitatively estimate the mobilization of trapped/slow moving oil mass in an idealized porous geometry, varying the Capillary number Ca and the wettability of the solid surfaces, expressed in terms of the contact angle $\vartheta$. The variation of these two parameters mimics the effect of the presence of surfactants and control both the adhesion to the walls and the degree of emulsification of the dispersed mass. High adhesion and low Ca lead to fewer free droplets, while their size seems to be weakly affected. Furthermore, the analysis shows that the adhesion highly influences the trapping process, inducing a thick continuous layer of slow mass stuck to the pore walls. However, departure from such configuration is attainable through even slight variations in the Ca, promoting re-mobilization of the adhered layers.

In conclusion, the LBM proves to be a suitable tool to tackle the complexity of the phenomena related to this kind of problems. Moreover, the analysis shows the capability of the proposed LBM model to give insights on the trapping/mobilization phenomenon. Therefore, the results open to a quantitative modelling of such complex processes in more complex and realistic geometries to further understand and improve the efficiency of processes like oil extraction and soil remediation.

Appendix A Additional information

Here some additional information about the model are presented. Nonetheless, the interested reader is still referred to Montessori et al. (2019) for an extended description of the adopted LBM model.

The three terms in Eq. 3 are calculated in the following way.

The first term is defined as:

$$\begin{aligned} (\Omega _a^k)^{(1)}=\omega _{eff}(f_a^{k,eq}-f_a^k) \end{aligned}$$

(A1)

It represents the standard single relaxation collisional operator (Succi 2001) with $\omega _{eff}=2/(6\overline{\nu }-1)$, called effective relaxation parameter, accounting for the effective kinematic viscosity of the two components $\overline{\nu }$, defined as in Eq. A2.

$$\begin{aligned} \frac{1}{\overline{\nu }}=\frac{\rho ^1}{(\rho ^1+\rho ^2)}\frac{1}{\nu ^1}+\frac{\rho ^2}{(\rho ^1+\rho ^2)}\frac{1}{\nu ^2} \end{aligned}$$

(A2)

The term $f_a^{k,eq}$ is the equilibrium PDF of the $k^{th}$ component and has the same form as in Eq. A3 if the two components have the same density ($\rho ^1/\rho ^2=1$).

$$\begin{aligned} f_{a}^{eq}=w_{a}\rho \left[ 1+\frac{(\textbf{e}_{a}\cdot \textbf{u})}{c_{s}^{2}}+\frac{(\textbf{e}_{a}\cdot \textbf{u})^{2}}{2c_{s}^{4}}-\frac{\textbf{u} \cdot \textbf{u}}{2c_{s}^{2}}\right] \end{aligned}$$

(A3)

being $w_a$ some weights (whose values depend on the velocity set chosen), $\textbf{u}$ the macroscopic velocity and $c_s$ the lattice sound speed.

The term $(\Omega _a^k)^{(2)}$ is the perturbation operator (Gunstensen et al. 1991), which is connected to the fluid surface tension $\sigma$; in particular, this operator must comply to the following constraints in order to reproduce the correct form of the stress tensor.

$$\begin{aligned}&\sum _{a=0}^M(\Omega _a^k)^{(2)}=0 \end{aligned}$$

(A4a)

$$\begin{aligned}&\sum _{k=0}^2\sum _{a=0}^M(\Omega _a^k)^{(2)}\textbf{e}_a=0 \end{aligned}$$

(A4b)

Furthermore, the perturbation operator is assumed having the shape (Leclaire et al. 2012):

$$\begin{aligned} (\Omega _a^k)^{(2)}=\frac{A^k}{2}|\nabla \psi |\Big [w_a\frac{(e_a\cdot \nabla \psi )^2}{|\nabla \psi |^2}-B_a\Big ] \end{aligned}$$

(A5)

Where $B_i$ is a set of parameters that must satisfy some isotropic constraints (Leclaire et al. 2017) and $A^k$ ($A^1=A^2$) is a parameter, which is connected to surface tension of the model $\sigma$ through the equation:

$$\begin{aligned} \sigma =\frac{4}{9}(A^1+A^2)\frac{1}{\omega _{eff}} \end{aligned}$$

(A6)

$(\Omega _a^k)^{(2)}$ generates an interfacial tension in compliance with the capillary-stress tensor of the Navier–Stokes equations for a multiphase system. However, it is not possible to ensure the immiscibility of different fluid using only such operator, hence the necessity of a further step, which is called recoloring step and is represented by the last term of Eq. 3.

$(\Omega _a^k)^{(3)}$ is the recolouring operator and minimizes the local diffusion between the two phases. For the two PDFs such operator takes the following form (Latva-Kokko and Rothman 2005a):

$$\begin{aligned} f_a^1&= \frac{\rho ^1}{\rho }f_a^*+\beta \frac{\rho ^1\rho ^2}{\rho ^2}\cos {\phi _a}f_a^{eq,0} \end{aligned}$$

(A7a)

$$\begin{aligned} f_a^2&= \frac{\rho ^2}{\rho }f_a^*-\beta \frac{\rho ^1\rho ^2}{\rho ^2}\cos {\phi _a}f_a^{eq,0} \end{aligned}$$

(A7b)

with $f_a^*=\sum _{k=1}^2f_a^{k,*}$ the set of post-perturbation PDFs, $\cos {\phi _a}$ the angle between the phase field gradient and the $a^{th}$ lattice direction and $f_a^{eq,0}=\sum _{k=1}^2f_a^k(\rho ,\textbf{u}=0)^{eq}$ the total equilibrium PDF calculated with zero velocity (Leclaire et al. 2012). The coefficient $\beta$ is a free parameter and controls the width of the interface (Latva-Kokko and Rothman 2005a).

It is worth specifying that in the adopted model the collision step (collision + perturbation + recoloring) are performed on two separate distribution as in Leclaire et al. (2012).

After the streaming step, a further regularization step, which consists in filtering out high-order non-hydrodynamic nodes that may emerge, is carried out (Mattila et al. 2017; Hegele Jr et al. 2018).

The distance d(x) is evaluated considering the phase field $\psi$ inside the domain: in the adopted model, a node is defined on the interface if $-0.05<\psi <0.05$. From an interface node, the algorithm searches if there is another one nearby along the surface normal $\textbf{n}$ in the distance $d_{max}$, which, for this class of models, is the typical thickness of the interface and does not scale with the resolution. The distance d(x) is then evaluated accordingly if such interface node is found.