1. Introduction
In chemical, geological and biological engineering applications, capillary trapping of non-aqueous phases (oil ganglia, bubbles or microorganisms) in porous media is a well-known phenomenon that commonly emerges in the aftermath of multiphase displacement (Macminn, Szulczewski & Juanes Reference Macminn, Szulczewski and Juanes2010; Lei et al. Reference Lei, Lu, Gong and Wang2023b
), component diffusion (Xu, Bonnecaze & Balhoff Reference Xu, Bonnecaze and Balhoff2017; Tan et al. Reference Tan, Banerjee, Shi, Tang, Abdel-Fattah and Squires Todd2021), or chemical reactions (Jiménez-Martínez et al. Reference Jiménez-Martínez, Hyman, Chen, Carey, Porter, Kang, Guthrie, George and Hari2020; Ladd & Szymczak Reference Ladd and Szymczak2021). In these porous systems, capillarity is hugely dominant over viscous, gravitational or inertial forces, where droplets and bubbles are confined and trapped in pores due to larger capillary force (Singh et al. Reference Singh, Jung, Brinkmann and Seemann2019; Lei et al. Reference Lei, Gong, Lu and Wang2024). Capillary trapping in porous media generally can be overcome by applying a larger viscous pressure drop (Datta et al. Reference Datta, Chiang, Ramakrishnan and Weitz2013). However, mobilization of trapped fluid from small, confined regions – such as dead-end pores of a porous medium – by pushing through a second fluid can be very difficult, if not impossible. In this paper, we maintain the dominance of capillary pressure in the dead-end pores of complex materials, but reverse the conventional capillary trapping via nanoparticle suspension to release the trapped fluid. Assuming that there exists some macroscopic surface modulation in the confined pore space, which is naturally occurring in engineering porous materials, such as rock and soil, we seed an invading liquid with nanoparticles that adsorb onto the surface spontaneously, creating a multiscale micro/nanostructure that allows the liquid film to penetrate the dead-end pore spontaneously even under weakly water-wet conditions. The established film on the multiscale surface results in a beneficial capillary pressure gradient to drive out the trapped liquid. The proposed strategy can also significantly enhance fluid displacement in a random porous medium. Our findings open a door for manipulating capillary pressure gradient direction to trap or release non-wetting fluid in complex geometry, and show a potential use in a wide range of applications, including microfluidic logic control (Cybulski, Garstecki & Grzybowski Reference Cybulski, Garstecki and Grzybowski2019), liquid-infused material fabrication (Yasuga et al. Reference Yasuga2021), cellular fluidics (Dudukovic et al. Reference Dudukovic, Fong, Gemeda, DeOtte, Cerón, Moran, Davis, Baker and Duoss2021), the remediation of contaminated soil and groundwater (Pak et al. Reference Pak, Luz, Tosco, Costa, Rosa and Archilha2020), and geological CO
$_2$
sequestration (Huppert & Neufeld Reference Huppert and Neufeld2014).
The fluid interfacial energy exerted on a droplet, bubble or particle by capillary pressure gradient is commonly reported to drive its motion in an open surface (Li et al. Reference Li, Qin, Shah, Ras Robin, Tian and Jokinen2016; Yadav, Hinch & Tirumkudulu Reference Yadav, Hinch and Tirumkudulu2019; Bintein et al. Reference Bintein, Lhuissier, Mongruel, Royon and Beysens2019; Wang et al. Reference Wang, Si, Dai, Li, Gao, Dong and Jiang2020), while few studies have shown that capillary pressure gradients can realize mobilization of trapped fluid in geometric-confined porous media. Combining connected water film with structure modulation in the confined region may create a favourable capillary pressure gradient to release the trapped fluid. However, super-hydrophilic conditions for film formation are hard to achieve and are harmful for compact displacement in porous media (Lei et al. Reference Lei, Lu, Liu and Wang2022, Reference Lei, Gong and Wang2023a ). Moreover, there is also a dilemma in combining the water film with the structure modulation in the dead-end pore. Thin liquid films on the structure modulation in the pore space can produce strong inside capillary pressure, while the corresponding convex solid surface of the macroscopic structure modulation will suppress the connected film formation because the liquid film is easier to break on the convex surface (Kovscek, Wong & Radke Reference Kovscek, Wong and Radke1993; Israelachvili Reference Israelachvili2011). Inspired by the multicurvature structures of the natural peristome surface of Nepenthes alata, which is adept at liquid harvesting and transport (Li et al. Reference Li, Yu, Zhou, Dong and Jiang2020; Hardt & McHale Reference Hardt and McHale2022), the integration of multiscale surfaces–combining nanoscale roughness with microscale structures–could potentially overcome the limitations associated with solid surface curvature. This multiscale surface introduces a significant nanoscale concave curvature on the microscale convex curvature, which mitigates the drawbacks of the underlying convex solid surface, and preserves the advantages of a concave solid surface.
Nanoparticles are ideal for creating nanostructures within complex porous geometries. Numerous experimental and theoretical studies have investigated the kinetics of nanoparticle suspensions at interfaces (Lei et al. Reference Lei, Lu and Wang2023c
; Nikolov, Wu & Wasan Reference Nikolov, Wu and Wasan2019). The adsorption of suspended nanoparticles on fluid–solid (Al-Anssari et al. Reference Al-Anssari, Barifcani, Wang, Maxim and Iglauer2016) or fluid–fluid interfaces (Isa et al. Reference Isa, Lucas, Wepf and Reimhult2011) may help to enhance capillarity by reducing interfacial tension (Binks Reference Binks2002), altering wettability (Al-Anssari et al. Reference Al-Anssari, Barifcani, Wang, Maxim and Iglauer2016), generating/stabilizing emulsions/foams (Sacanna, Kegel & Philipse Reference Sacanna, Kegel and Philipse2007), and reacting with chemical components (Xu, Agrawal & Darugar Reference Xu, Agrawal and Darugar2018; Pak et al. Reference Pak, Luz, Tosco, Costa, Rosa and Archilha2020). A common effect of nanoparticles in a three-phase contact region is the formation of suspended nanoparticle layers in the confined three-phase edge region, which increases the disjoining pressure, driving the inward movement of the three-phase interface (Wasan & Nikolov Reference Wasan and Nikolov2003). However, this disjoining pressure effect is significant only under super-hydrophilic conditions (water contact angle
$\lt5^\circ$
) and ultra-low interfacial tension (
$\lt$
1 mN m–1), and the volume fraction of particles must be greater than 10 vol.
$\%$
(Wasan & Nikolov Reference Wasan and Nikolov2003; Nikolov et al. Reference Nikolov, Wu and Wasan2019). The applicability of this unique effect in enhancing multiphase displacement in porous media has been discussed in our previous review paper (Lei et al. Reference Lei, Lu and Wang2023c
). Furthermore, these effects typically overlook the influence of complex geometric confinement, a common feature in various porous systems, such as rocks, soils, packed beds and membranes, where pore surface structure, pore angularity and pore connectivity make interfacial phenomena more complex.
Through microfluidic experiments and numerical simulations, we demonstrate that nanoparticle suspensions can initiate a robust water film and reverse stable capillary trapping, thereby releasing non-aqueous fluid from dead-end pores. Our findings reveal a multiscale capillarity between nanoparticle-induced nanoscopic structures and the inherent microscopic dead-end formations. This multiscale surface allows for the spontaneous generation of a stable hydrodynamic film via capillary condensation within the nanoscopic structures, and maintains the connectivity of the hydrodynamic film. Utilizing this connected hydrodynamic film as a prerequisite, trapped non-aqueous fluid in specific dead-end geometries can be completely released by the capillary force gradient. The dynamics and evolution of this self-removal process, driven by a capillary pressure gradient, can be predicted accurately using mathematical models that align well with experimental data. Additionally, a regime map based on pore angularity and microscopic roughness ratio can forecast whether trapped fluid can be entirely released. This reversing capillary trapping phenomenon is general in multiphase flow in porous media, which holds significant implications for enhancing displacement in random porous media, potentially leading to novel advancements in controlling multiphase flow across various engineering applications.
2. Methods
2.1. Experimental set-up and materials
Figure 1(a) illustrates the microfluidic platform designed to investigate pore-scale interfacial phenomena in dead-end microchannels, and macroscopic invasion patterns in microfluidic porous media. The cleaned microchip was assembled in a holder and placed under a fluorescence microscope (Nikon SMZ18) equipped with a high-speed camera (Nikon DS-Ri2). The injection system comprises a microfluidic pressure controller (Fluigent MFCS-8C), a flow rate sensor (Fluigent FRP), a pressure transducer (Validyne P61) and two reservoirs. The displacement processes were visualized and quantified by the in-house image processing and data analysis software. For the microfluidic dead-end experiments, spatial resolution 0.96
$\unicode{x03BC}\text{m}\ \text{pixel}^{-1}$
and time resolution 1 min were used. In the porous media displacement experiments, spatial resolution 1.92
$\unicode{x03BC}\text{m}\ \text{pixel}^{-1}$
and time resolution 1 s were selected to ensure quasi-static conditions, allowing analysis of interfacial phenomena at the microscale and invasion patterns at the centimetre scale.
The microfluidic platform and microfluidic chips. (a) A schematic diagram of the experimental system. (b) Glass microfluidic chips for dead-end pores and silicon microfluidic chips for porous media.
To observe the trapping of non-aqueous fluid in dead-end pores and the evolution of fluid–fluid interfaces, we conducted microfluidic experiments using structured microchannels containing several dead-end pores. Each dead-end structure was connected to a main flow channel for pore-scale studies. To obtain representative displacement results, the microfluidic porous media were generated using the quartet structure generation set method, based on statistical information about rock obtained from a real geological system, as introduced in our previous work (Lei et al. Reference Lei, Liu, Xie, Yang, Wu and Wang2020). All dead-end structures were fabricated from glass, while the porous structures were made from silicon due to the small and dense geometries of the microfluidic chips used in this study (figure 1
b). The structures were created on the substrates using microfabrication techniques, including standard photolithography, HF acid etching for glass substrates, and inductively coupled plasma-deep reactive ion etching (ICP-DRIE) for silicon wafers (Lei et al. Reference Lei, Lu and Wang2023c
). The microchips were sealed by anodically bonding a thin glass layer to ensure proper visualization. To preserve the nanoparticle adsorption effect, the etched silicon wafer was heated to 1000
$^\circ$
C under oxygen to create a uniform 100 nm silicon dioxide layer on the etched porous structures. The flow rate was controlled using a feedback-controlled pump (MFCS-EZ, Fluigent) equipped with a precise flow meter (Flowunit, Fluigent).
The non-aqueous fluid was decane (95
$\%$
, Macklin), dyed with Sudan III (100 ppm, Macklin) for dead-end experiments, or Nile Red (100 ppm, Macklin) for porous media displacement experiments to enhance visualization. A colloidal dispersion of silicon oxide in water (liquid–liquid interfacial tension
${\sigma } = 44.3$
mN m–1, nominal diameter 20 nm, salinity 1.350 g l–1, and pH 9.7) was used as the aqueous phase. An aqueous surfactant solution (0.01 wt.
$\%$
sodium dodecyl sulphate (SDS) with
$\sigma=43.1$
mN m–1, salinity 1.350 g l
$^{-3}$
by NaCl, and pH 9.3) was used as a comparison group to minimize the influence of interfacial tension, salinity and pH. The colloidal dispersion of silicon (IV) oxide in water (40 wt.
$\%$
, Alfa Aesar) was diluted with deionized water to prepare a suspension containing 10 vol.
$\%$
silica nanoparticles. The aqueous surfactant solution (SDS, VWR) was diluted in brine water to 0.01 wt.
$\%$
, with brine water prepared using sodium chloride (NaCl) and sodium hydroxide (NaOH) to match the salinity and pH. The system exhibited weakly water-wet wettability, with contact angle
$\theta = {48}^{\circ }$
, where corner flow could be ignored. Corner flows typically arise when
${\theta } \lt {90}^{\circ } - {\varphi }/2$
(Concus & Finn Reference Concus and Finn1969), due to geometrical features such as corners and wedges, where
$\varphi$
is the corner angle of the cross-section in a microchannel (here,
$\varphi = {90}^{\circ }$
).
2.2. Image analysis
Ganglia size distributions of the non-aqueous fluid and surface curvatures were analysed using the image analysis software ImageJ. Since all fluid ganglia share the same curvature in the dimension perpendicular to the microscopic image plane, we calculated the differential capillary pressures based only on the curvature in the horizontal dimension (the plane of the microscopic images). It is worth noting that while ImageJ can automatically calculate curvature for large droplets, the curvatures of all surfaces (including angular areas and throats) were measured manually five times independently and then averaged during the early stage of self-removal. The MATLAB Image Processing Toolbox, along with a custom-developed program, was used to analyse the displacement in porous media. Both the saturation and the fluid topology are quantified via the Euler characteristics
$\chi ={\beta _0}-{\beta _1}+{\beta _2}$
, where
$\beta _0$
is the number of objects,
$\beta _1$
is the number of redundant loops, and
$\beta _2$
is the number of cavities ignored in experiments. The shape factor or circularity is
$c={4\unicode{x03C0} A}/{p^2}$
(where
$A$
is the ganglia area, and
$p$
is the ganglia perimeter). The morphology type of the residual oil was classified as droplets (
$\chi =1$
and
$0.5\lt c\leqslant 1$
), columns (
$\chi =1$
and
$0\lt c\leqslant 0.5$
) and cluster (
$\chi \lt 1$
) (Lei et al. Reference Lei, Lu, Liu and Wang2022).
2.3. Numerical solution of the augmented Young–Laplace equation
The finite-difference numerical solution of the augmented Young–Laplace equation was obtained using successive over-relaxation (SOR) for the water film configuration. The augmented Young–Laplace equation is given by (Morrow Reference Morrow1990)
\begin{equation} {P_c} = \sigma \kappa + \Pi \left ( h \right ) = \sigma \left ( {\frac {1}{{{R_1}}} + \frac {1}{{{R_2}}}} \right ) + \Pi \left ( h \right ), \end{equation}
where
$R_1$
and
$R_2$
are the principal radii of curvature at each point on the fluid–fluid interface
$z=w(x,y)$
. The term
$\Pi ( h )$
represents the disjoining pressure of the thin film, which depends on the film thickness
$h$
and can be determined through extended DLVO interaction force measurements and calculations, as detailed in § 3.2.
The local mean curvature of the menisci geometries,
$\kappa$
, can be derived from the derivatives for the coordinate directions
$x$
and
$y$
in Cartesian coordinates of the rough structure plane (Kibbey Reference Kibbey2013):
\begin{equation} \kappa = \frac {1}{{{R_1}}} + \frac {1}{{{R_2}}} = \frac {{\left ( {1 + {{\left ( {\frac {{\partial w}}{{\partial x}}} \right )}^2}} \right )\frac {{{\partial ^2}w}}{{\partial {y^2}}} - 2\frac {{\partial w}}{{\partial x}}\frac {{\partial w}}{{\partial y}}\frac {{{\partial ^2}w}}{{\partial x\,\partial y}} + \left ( {1 + {{\left ( {\frac {{\partial w}}{{\partial y}}} \right )}^2}} \right )\frac {{{\partial ^2}w}}{{\partial {x^2}}}}}{{{{\left ( {1 + {{\left ( {\frac {{\partial w}}{{\partial x}}} \right )}^2} + {{\left ( {\frac {{\partial w}}{{\partial y}}} \right )}^2}} \right )}^{{3 / 2}}}}}. \end{equation}
It is important to note that while the choice of the coordinate system does not influence the results, it is essential to ensure that
$z = w(x, y)$
can be uniquely determined by
$x$
and
$y$
. This condition is also consistent with the characteristics of the water film.
By substituting (2.2) into (2.1), we solved the augmented Young–Laplace equation using the finite-difference method with SOR. Neumann boundary conditions were applied at the fluid–fluid interface edge, where the boundary slope of the calculated fluid–fluid interface was set equal to the nearby slope of the fluid–fluid interface. After the initial fluid–fluid interface was established, the actual interface was iteratively solved for all fluid–fluid interface points not in contact with the solid surface, using the equation
\begin{equation} \frac {{\partial w}}{{\partial \phi }} = \kappa - {{{P_c}} / \sigma }, \end{equation}
where
$\phi$
is an iteration parameter (with units of mesh size squared). Each step was solved using SOR. As
$\kappa$
approaches
$P_c$
on the water interface, the left-hand side of (2.3) approaches zero, resulting in an interface that satisfies the augmented Young–Laplace equation (2.1). The accuracy of the numerical algorithm is benchmarked and validated in Appendix A.
2.4. Multiphase lattice Boltzmann simulations
The experiments and simulations flow system is under normal temperature and pressure conditions. The governing equations for the isothermal incompressible two-phase immiscible Newtonian fluids are
\begin{equation} \frac {{\partial \left ( {\rho \boldsymbol {u}} \right )}}{{\partial t}} + \boldsymbol{\nabla} \boldsymbol{\cdot} \left ( {\rho \boldsymbol {uu}} \right ) = - \boldsymbol{\nabla} p + \boldsymbol{\nabla} \boldsymbol{\cdot} \left ( {2\mu \boldsymbol {E}} \right ) + {\boldsymbol {f}_s}, \end{equation}
\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol {u} = 0, \end{equation}
where
$t$
designates time,
$\rho$
is the density,
$\mu$
is the dynamic viscosity,
$\boldsymbol {E} = ({1}/{2})$
$({\boldsymbol{\nabla} \boldsymbol {u} + \boldsymbol{\nabla} ^{\rm T}\boldsymbol {u}})$
is the rate-of-strain tensor,
$\boldsymbol u$
is the vector of fluid velocity,
$p$
is the pressure, and
$\boldsymbol {f}_s$
is the volumetric surface tension force that describes the effect of the Laplace pressure.
Based on mesoscopic kinetic equations, the lattice Boltzmann method (LBM) is particularly effective for solving the above equations with complex boundary conditions and multiphase interfaces (Krüger et al. Reference Krüger, Kusumaatmaja, Kuzmin, Shardt, Silva and Viggen2017; Xie, Lei & Wang Reference Xie, Lei and Wang2018; Xie & Balhoff Reference Xie and Balhoff2021; Zheng et al. Reference Zheng, Ju and Wang2018, Reference Zheng, Lei, Ju and Wang2021). There are four basic multiphase LBMs, which include the colour gradient model, the pseudo-potential lattice Boltzmann method (SC-LBM), the free-energy model, and the phase-field model (Aidun & Clausen Reference Aidun and Clausen2010). Among them, the SC-LBM is one of the most popular multiphase flow methods due to its easy implementation. In this study, we employed SC-LBM to simulate the self-removal process of nonaqueous fluid from the dead end. The general framework of the model follows our previous work (Zheng et al. Reference Zheng, Ju and Wang2018, Reference Zheng, Lei, Ju and Wang2021), utilizing the D3Q19 square lattice with explicit forcing terms and the multiple relaxation time method. To meet the computational demands of simulating a three-dimensional dead-end structure with multiscale micro/nanostructures, we developed an in-house CUDA-LBM code (Yang et al. Reference Yang, Chen, Chen and Wang2023), which achieved 803 million lattice updates per second (MLUPS) on a single NVIDIA V100 GPU, providing a speed-up of approximately 2000 times over serial implementation. Despite this efficiency, simulating a single case of fluid–fluid interface evolution within the dead-end pore still required over six months on an NVIDIA V100 GPU. To mimic the nanoscopic structure induced by nanoparticle adsorption, pillar-like structures were generated on the surface of the dead-end pore. Details of the simulation method and benchmark validation are provided in Appendix B.
Releasing non-aqueous fluid from a dead end via a nanoparticle suspension under weakly water-wet conditions. (a) Schematic showing the typical dead-end structure of the microchip. (b) Experimentally observed spontaneous release of trapped non-aqueous fluid from the dead-end structure, and corresponding colour-modified images. (c) Representative results of fluid release from the various dead-end systems via nanoparticle suspension and comparison with aqueous surfactant solutions at the final stage (at 4 h).
Multiscale surface, disjoining pressure, curvature radius-related film thickness, and the transition between hydrodynamic film and adsorption film. (a) Scanning electron microscope images of the dead-end structure and the microscopic convex surface. (b) Atomic force microscopy images of nanoscopic structures formed by nanoparticle adsorption and original surface without nanoparticle adsorption. (c) Disjoining pressure isotherm for wetting films on solids under our experimental conditions. (d) Schematic diagram of hydrodynamic film rupture from hydrodynamic film to adsorption film due to a small curvature radius of the convex surface, and capillary condensation from adsorption film to hydrodynamic film due to a small curvature radius of the concave surface.
3. Results and discussion
3.1. Pore-scale dynamics of spontaneous fluid removal
To impose non-aqueous fluid trapped in a dead-end pore, we used a microfluidic approach, in which the entrance of the dead-end structure was connected to a main flow channel, and the typical dead-end structure contained a regular polygon pore space with a microscopic rough structure (figure 2
a). The initial capillary pressure
${P}_{c}$
in the dead end was controlled as
${{{P}}_c} = 2\sigma ( { ({1}/{W}) + ({1}/{D}})) = 6.2$
kPa in all the dead-end systems, where the throat width
$W$
was 50
$\unicode{x03BC}\rm m$
, and the depth
$D$
was 20
$\unicode{x03BC}\rm m$
. During these microfluidic experiments, we surprisingly observed the spontaneous release of non-aqueous fluid from the dead end via a nanoparticle suspension, as shown in figure 2(b). To clarify the evolution and phase distribution of the fluid–fluid interface, the water phase was reassigned to light blue, the oil phase to orange, and the solid matrix to grey. Compared to the conventional capillary trapping via the aqueous surfactant fluid, the nanoparticle suspension can reverse this capillary trapping and transport into the dead end as a film state to release non-aqueous fluid from various dead-end geometries, as shown in figure 2(c). Notably, the trapped non-aqueous fluid can be entirely released from some dead-end systems, which are controlled by the pore shape and the microscopic structure. Appendix C and supplementary movies 1–5 available at https://doi.org/10.1017/jfm.2025.53 present more detailed experimental results. These results agree with our hypothesis that a nanoparticle-stimulated water film combined with the dead-end structure drives trapped non-aqueous fluid out from dead-end pores, leading to this counterintuitive capillary phenomenon.
3.2. Underlying microscopic physics of nanoparticle-enhanced film flow
We analysed the initial water film formation and configuration on the dead-end structure to elucidate the microscopic mechanism of this anomalous capillary phenomenon. The microstructure in the dead-end structure and the nanoscopic roughness formed by nanoparticle adsorption were characterized by scanning electron microscope and atomic force microscopy (AFM), respectively. The curvature radius of the microscopic convex surface in the dead end was approximately 100 nm (figure 3
a). Quantitative roughness characterizations of the original surface (without nanoparticle adsorption) and the nanoparticle-adsorbed surface, as measured by AFM, are presented in table 1. The in situ nanoscopic structure formed by nanoparticle adsorption on the glass substrate was characterized as the root mean square roughness
$R_{rms}$
, which is 6.8 nm, and the peak-to-valley roughness
$R_t$
is 72.5 nm (figure 3
b).
The curvature distribution in micro/nanostructures (figures 3
a and 3
b) and disjoining pressure isotherm (figure 3
c) will determine film formation and evolution. We characterize the extended DLVO interaction forces and disjoining pressure under the experimental environment. Interfacial forces acting between the solid and aqueous fluid surfaces present in a non-aqueous fluid were estimated using the extended DLVO theory (Morrow Reference Morrow1990; Israelachvili Reference Israelachvili2011). Free energy
$W$
is the summation of van der Waals (VDW), electrostatic or electrical double layer (EDL) and structural (S) free energies as a function of film thickness
$h$
:
\begin{equation} W ( h ) = {W_{EDL}} ( h ) + {W_{VDW}}( h ) + {W_S} ( h ). \end{equation}
The derivative of these interaction energies per unit area concerning the separation distance (film thickness
$h$
) in the direction normal to that of the interacting surfaces is known as the disjoining pressure (
$\Pi$
):
\begin{equation} \Pi = - \frac{{\partial ( {{W_{EDL}} ( h ) + {W_{VDW}} ( h ) + {W_S} ( h )} )}}{{\partial h}}. \end{equation}
The electrostatic force, called the EDL force, is caused by the
$\zeta$
-potential at the surfaces. The surface of solid and non-aqueous fluid will often be charged when in contact with the aqueous fluid due to the ionization and dissociation of surface groups on the solid and non-aqueous fluid surface. The EDL free energy can be obtained by solving the Poisson–Boltzmann equation (PBE). Based on the analytical solution of the linearized PBE under the constant potential boundary condition, we can derive the EDL interaction energy as (Israelachvili Reference Israelachvili2011; Lei et al. Reference Lei, Lu and Wang2023c
)
\begin{equation} {W_{EDL}} ( h ) = \frac {{\varepsilon {\varepsilon _a}k\left [ {2{\zeta _1}{\zeta _2} - \left ( {\zeta _1^2 + \zeta _2^2} \right ){{e}^{( { - kh} )}}} \right ]}}{{2\sinh ( {kh} )}}, \end{equation}
where
$\zeta _1$
and
$\zeta _2$
are the
$\zeta$
-potential on the solid–aqueous fluid and non-aqueous–aqueous fluid interfaces, respectively. Also,
${k^{ - 1}} = \sqrt { ( {\varepsilon {\varepsilon _a}{K_B}T} ) / ( {2{N_A}{{e}^2}I} )}$
is the Debye length, wherein
$I=0.5\sum {z_i^2\rho _i}$
is the ionic strength,
$\rho _i$
is the ion concentration (mol kg
$^{-1}$
), and
$z_i$
is the number of charges carried by the ion.
The van der Waals force dominates when the film thickness (
$h$
) is small. The van der Waals energy (
$W_{VDW}$
) is usually present and originates due to the dipole-dipole, dipole-induced-diple and dispersion forces (Morrow Reference Morrow1990; Lei et al. Reference Lei, Lu and Wang2023c
):
\begin{equation} W_{VDW} (h )=-{A}/{ (12\unicode{x03C0} h^2 )\ }, \end{equation}
where
$A = ({3}/{4}){K_B}T ( ({ {{{\varepsilon _n} - {\varepsilon _a}}})/({{{\varepsilon _n} + {\varepsilon _a}}}}) ) ( ({ {{{\varepsilon _s} - {\varepsilon _a}}})/({{{\varepsilon _s} + {\varepsilon _a}}}}))$
is the Hamaker constant as a function of the dielectric permittivities of different phases–solid (
$\varepsilon _s$
), aqueous fluid (
$\varepsilon _a$
) and non-aqueous fluid (
$\varepsilon _n$
)–which varies mainly with temperature. Here,
$K_B$
is the Boltzmann constant, and
$T$
is the Kelvin temperature.
Roughness characterization.
The structural force is an interaction between two surfaces when the separation or film thickness is very small. It is often called the solvation force when the medium is a solvent, and the hydration force if the solvent is water. The structural interaction is calculated from the expression (Morrow Reference Morrow1990; Lei et al. Reference Lei, Lu and Wang2023c )
\begin{equation} W_S (h)=h_sA_s\,{\rm e}^{-{h}/{h_s}}, \end{equation}
where
$A_s$
is the empirical coefficient, and
$h_s$
is the characteristic decay length for the exponential model.
All the parameters for the extended DLVO forces calculation in this study are presented in table 2. According to the DLVO curves in this study (figure 3
c), the adsorption film is defined as an ultra-thin film (
$\lt$
0.56 nm) with several water molecules dominated by structural force. In comparison, the hydrodynamic film is defined as a thick film (
$\gt$
2.0 nm) dominated by van der Waals forces and EDL forces. It is worth noting that a uniform molecular-level adsorption film is preset on the wall based on DLVO theory for hydrophilic surfaces, but the subsequent fluid transport for self-removal processes can only rely on hydrodynamic film.
Parameters for the extended DLVO model.
[a] The
$\zeta$
-potential of the silica particles and aqueous fluid was measured using Zetasizer Nano-ZS (Malvern Instruments). A colloidal dispersion of silicon (IV) oxide in water (40 wt.%, Alfa Aesar) was diluted with deionized water to prepare the 10 vol.% silica nanoparticle suspension used as the measuring solution.
[b] To improve the accuracy of the measurements and ensure higher purity of the system, we first produced finely dispersed, evenly sized decane droplets in the pure water, following the methods proposed by Marinova et al. (Reference Marinova, Alargova, Denkov, Velev, Petsev, Ivanov and Borwankar1996). Then the solution conditions, such as salinity and pH, are adjusted to be the same as in the nanoparticle suspension; the dispersion did not separate back for a few tens of minutes. Finally, Zetasizer Nano-ZS (Malvern Instruments) was used to measure the non-aqueous/aqueous fluid
$\zeta$
-potential.
[c] The ion concentration of the nanoparticle suspension is measured by inductively coupled plasma mass spectrometry (ICP-MS, ThermoFisher), with the content of each element being Na 0.02320 mol l–1, S 0.001533 mol l–1, Al 0.0008107 mol l–1, and Cl 0.0295926 mol l–1.
[d] The inverse of the Debye length
$k$
is calculated based on
${k^{ - 1}} = \sqrt { ( {\varepsilon {\varepsilon _a}{K_B}T} ) / ( {2{N_A}{e^2}I} )}$
; all the parameters can be acquired from this table.
The conversion from the molecular-level adsorption film to hydrodynamic film is difficult due to the limitation of disjoining pressures
$\Pi ^{max}$
and
$\Pi ^{min}$
(figure 3
c). Multiscale surfaces by generating nanoscopic structures on the microscopic convex surfaces may promote this conversion (figure 3
d). When the adsorption film is on the concave surface and the minimum disjoining pressure cannot balance the capillary pressure
$P_c\lt \sigma \kappa +\Pi ^{min}$
because of
${r_{cc}} \lt {r_{cc\_cri}}$
, where the minimum disjoining pressure is
$\Pi ^{min}=-430$
kPa, the critical value of concave curvature radius is
${r_{cc\_cri}} = {{2\sigma } / { ( {{P_c} - {\Pi ^{min }}} )}}= 200$
nm, and
$\kappa$
is the surface curvature, it can trigger the formation of a hydrodynamic film on the nanoscopic concave surface from the stable molecular-level adsorption film due to capillary condensation (Tuller, Or & Dudley Reference Tuller, Or and Dudley1999; Anand et al. Reference Anand, Ghosh, Aabdin, Koneti, Xu, Holsteyns and Mirsaidov2021). When the hydrodynamic film is on the smooth convex curvature and the maximum disjoining pressure cannot balance the capillary pressure
$P_c\gt \sigma \kappa +\Pi ^{max}$
because of
${r_{cv}} \lt {r_{cv\_cri}}$
, where the maximum disjoining pressure is
$\Pi ^{max}= 63.5$
kPa and the critical value of convex curvature radius is
${r_{cv\_cri}} = {{ - 2\sigma } / { ( {{P_c} - {\Pi ^{max}}} )}}$
= 1.55
$\unicode{x03BC}$
m, there only exists stable molecular-level adsorption film due to hydrodynamic film rupture; see figure 3(d).
In our microfluidics experiments, on the microscopic smooth convex surface, there only exists a stable molecular-level adsorption film because
${r_{cv}} = 100\ {\rm nm}<{r_{cv\_cir}} = 1.55$
$\unicode{x03BC}\rm m$
(figures 3
a,d). However, nanoparticles can adsorb on the above microscopic convex surface and reverse this stable state because the concave curvature radius of nanoscopic structure
${r_{cc}}\lt {r_{cc\_cri}}= 200$
nm (figures 3
b,d); thus the adsorption film is transformed into a hydrodynamic film, and the robust connectivity of the hydrodynamic water film can be maintained by nanoscopic structures rather than destroyed by the underlying microscopic convex surface.
Nanoparticle-stimulated film formation and development on the multiscale surface. (a) The calculated aqueous–non-aqueous fluid interface on the multiscale micro/nanostructures under different capillary pressures; the blue colour is the hydrodynamic water film. (b) The calculated average thickness of the water film as a function of capillary pressure for the smooth microscopic surfaces and the multiscale micro/nanosurfaces with various curvature radii (
${r}_{cv}$
).
To intuitively show the formation of hydrodynamic water film on this multiscale micro/nanostructure, we further solved the augmented Young–Laplace equation using the finite difference method with SOR by setting the molecular-level adsorption film on the smooth and multiscale surfaces, respectively. Figure 4(a) compares the hydrodynamic water film configuration on the smooth microscopic convex surface and the multiscale micro/nanosurface under various capillary pressure conditions, where the multiscale micro/nanosurface is designed by adding nanoparticle adsorption morphology on the smooth microscopic convex surface, as shown in figures 3(a) and 3(b). On smooth microscopic convex surfaces, there only exists a molecular-level adsorption film regardless of the capillary pressure conditions, therefore the hydrodynamic film is disconnected or does not exist. However, on the multiscale micro/nanosurface, even at a high capillary pressure
$P_c = 62$
kPa (10 times the capillary pressure of current experiments,
${{{P}}_c} = 2\sigma ( { ({1}/{W}) + ({1}/{D}})) = 6.2$
kPa), the connected hydrodynamic water film will form spontaneously from the molecular-level adsorption film, and remain robust on the convex surface. Nanoscopic structures facilitate the transition from adsorption film to hydrodynamic film even on a macroscopic convex solid surface, which is impossible for a smooth convex solid surface (figure 4
a). The hydrodynamic film on the multiscale micro/nanostructure is orders of magnitude thicker than the adsorption film on the convex surface of the smooth microstructure, as shown in figure 4(b).
With the robust hydrodynamic film as a prerequisite, we analysed the subsequent release process from the dead end by direct numerical simulations via SC-LBM. The square dead-end system is selected as a representative structure for simulations. The simulations in dead-end structures with or without nanoscopic structures are compared. The dead-end system is approximately 500
$\unicode{x03BC}\rm m$
(figure 2
a). The mesh size of the dead-end system is 5 nm, with typical mesh count approximately
$8\times 10^{13}$
. The system dimensions are
${10}^5$
length,
$2\times {10}^5$
width, and
$4\times 10^{3}$
height, significantly exceeding current computational capabilities. As a result, the structure is scaled and simplified (figure 5
a). When the thickness of the structures is 30 lattices, the mesh size for each case is approximately
$3\times {10}^7$
. Although the mesh size is reduced, the time step needed is still tremendous. Over
$2\times {10}^8$
steps are needed to study the self-removal process. To our knowledge, no serial LBM implementation can support such mesh size and time step requirements. To accelerate LBM simulation, the LBM code is implemented via CUDA, which achieved 803 MLUPs on the NVIDIA V100 GPU. The stability of the hydrodynamic water film was compared in the dead end without and with the nanoparticle-modified surface (figures 5
b and 5
c). The stable hydrodynamic water film can only be maintained on the dead-end surface with nanoscopic structures. The interfacial curvature difference induces a capillary pressure gradient, gradually driving the trapped oil out (figure 5
d). Figure 5(d) shows the fluid–fluid interface evolution during the self-removal process. Supplementary movie 6 presents the simulated self-removal process.
3.3. Mathematical model and self-removal phase diagram
After the hydrodynamic film was formed, a mathematical model was developed to understand this anomalous capillary phenomenon better, and predict it. To develop a simplified but useful model, we assumed the following. (i) Nanoparticles can trigger a hydrodynamic water film and maintain its connectivity, which has been demonstrated in the above microscopic physics analysis (§ 3.2). (ii) The spontaneous release process is quasi-static with negligible viscous force, inertial force and gravity, because the time scale of the self-removal process is approximately 4 h, the capillary number is
${{Ca}} = {{\mu U} / \sigma } \sim {10^{ - 7}}$
, and the Bond number is
$Bd = {{\rho g{r^2}} / \sigma } \ll 1$
, where
$\mu$
is the fluid viscosity,
$U$
is the characteristic flow velocity,
$\rho$
is the fluid density,
$g$
is the gravitational acceleration, and
$r$
is the characteristic pore size. (iii) The wettability of microchips is weakly water-wet without corner flow because
$\theta = {48^ \circ } \gt {45^ \circ }$
. (iv) The roughness of the surface induced by the micro-manufacture is sufficiently small (
$R_{rms}=1.5$
nm,
$R_t=8.7$
nm); see table 1. Based on the above assumptions, the mathematical expression for the self-removal process of trapped non-aqueous fluid from dead-end structures can be derived based on mass conservation and capillary pressure gradient direction. Defining an appropriate indicator to characterize the fluid mobilization direction is helpful. Here, we introduce
${R^*} = {{ ( {{R_{out}} - {R_{in}}} )} / {{R_t}}}$
, the dimensionless difference between the inner and outer curvature radii, as the index of capillary pressure gradient direction, where
$R_{in}$
,
$R_{out}$
and
$R_t$
are the minimum interfacial curvature radius inside the dead end, the interfacial curvature radius outside the dead end, and the curvature radius of the released non-aqueous fluid with dead-end area
$A$
in the main channel, respectively (figure 6
a). Here,
${R^*} \gt 0$
means that the release process continues to be driven by capillary pressure gradient until the final complete release state
${R^*}=1$
, and
${R^*} \leqslant 0$
means that the release process stops due to capillary trapping, as shown in figure 6(b).
Numerical simulation results of hydrodynamic film connectivity and fluid–fluid evolution during self-removal. (a) Geometrical information about the dead-end system in the simulation, and local view of the nanoscopic structure generated on the surface (top) compared with a smooth surface (bottom). (b,c) Comparison of direct numerical simulation results of the hydrodynamic film connectivity in a dead end (b) without and (c) with nanoscopic structures on the microscopic surface. Inset: film distribution in the local region; the non-aqueous fluid is empty to present the film distribution more clearly. (d) Direct numerical simulation results for the fluid–fluid interface evolution in the dead end.
Physical models of dead-end structures. (a) Structure information and corresponding fluid–fluid interface evolution stages in typical dead-end structures. (b) The trapping or release state at the final stage.
As presented in § 3.1, the self-removal process goes through three stages, classified as the film-driven stage, the angular-radius-driven stage and the droplet-radius-driven stage (figure 6
a). According to the mass balance, the aqueous fluid saturation
$S_a$
as dimensionless time in these three stages can be written as
\begin{equation} S_a=\frac {\unicode{x03C0} R_{out}^2}{A}. \end{equation}
(i) In the film-driven stage (
$0 \leqslant {S_a} \lt {A^ * }$
), the film grows in the microscopic rough structure. The dimensionless difference between the inner and outer curvature radii as the capillary pressure difference
$R^\ast$
can be described as
\begin{equation} {R^ * } = \frac {{{R_{out}} - {R_{in}}}}{{{R_t}}} \approx \frac {{{R_{out}}}}{{\sqrt {{A / \unicode{x03C0} } } }} = \sqrt {{S_a}}, \end{equation}
where
${A^ * } = {{{A_r}} / A}$
is microscopic rough area factor, and
$A_r$
is the microscopic rough area of the dead-end structures. Here,
$R_{in}\rightarrow 0$
compared with
$R_{out}$
in the film-driven stage.
(ii) In the angular-radius-driven stage (
${A^ * } \leqslant {S_a} \lt {A^ * } + ({{{F_n}}}/{{2\unicode{x03C0} })}$
), when the water film detaches from the microscopic rough area, the angularity of the dead-end will dominate the fluid removal process. Here, we assume that the development of each corner of the polygon is similar, and then we can derive
$R_{in}$
in the angular-radius-driven stage:
\begin{equation} {R_{in}} = \sqrt {\frac {\unicode{x03C0} }{{ {n\tan ( {{\unicode{x03C0} / n}} ) - \unicode{x03C0} } }}\,\frac {{\left ( {{S_a} - {A^ * }} \right )A\;}}{\unicode{x03C0} }}, \end{equation}
where
$n$
is the number of regular polygon sides.
Thus we derive
$R^\ast$
as
\begin{equation} {R^*} = \frac {{{R_{out}} - {R_{in}}}}{{{R_t}}} = \sqrt {{S_a}} \left ( {1 - \sqrt {\frac {{\unicode{x03C0} ( {{S_a} - {A^*}} )}}{{{S_a}\left ( {n\tan ( {{\unicode{x03C0} / n}} ) - \unicode{x03C0} } \right )}}} } \right ) = \sqrt {{S_a}} - \sqrt {\frac {{\unicode{x03C0} \left ( {{S_a} - {A^*}} \right )}}{{n\tan ( {{\unicode{x03C0} / n}} ) - \unicode{x03C0} }}}. \end{equation}
Considering the definition of the pore-shape-dependent angularity factor
$F_n$
(Tuller et al. Reference Tuller, Or and Dudley1999), we have
\begin{equation} {F_n} = n\tan ( {{\unicode{x03C0} / n}} ) - \unicode{x03C0} = \sum_{i = 1}^n \left ( \frac {1}{\tan ( {{\alpha / 2}} )} - \frac {\unicode{x03C0} - \alpha }{2}\right), \end{equation}
where
$\alpha$
is the corner angle of the
$n$
-regular polygonal pore, with
$\alpha=\unicode{x03C0}$
for the circular pore (
$n=\infty$
,
$F_n=0$
), and
$\alpha={\unicode{x03C0} }/3$
for the triangular pore (
$n=3$
,
$F_n=2.1$
) .
Combining (3.9) and (3.10),
$R^\ast$
in the angular-radius-driven stage can be described as
\begin{equation} {R^ * } = \sqrt {{S_a}} - \sqrt {{\unicode{x03C0} / {{F_n}}} ( {{S_a} - {A^ * }} )}. \end{equation}
(iii) In the droplet-radius-driven stage (
${A^ * } + ({{{F_n}}}/{{2\unicode{x03C0} }}) \lt {S_a} \leqslant 1$
), when the non-aqueous fluid forms a circular shape within the dead-end structure without interacting with the microscopic roughness or angularity,
$R^\ast$
can be described as the radius difference inside and outside the dead-end:
\begin{equation} {R^ * } = \frac {{{R_{out}} - {R_{in}}}}{{{R_t}}} = \frac {{{R_{out}} - {R_{in}}}}{{\sqrt {{A / \unicode{x03C0} }} }} = \frac {{{R_{out}} - {R_{in}}}}{{{{{R_{out}}} / {\sqrt {{S_a}} }}}} = \left ( {1 - \frac {{{R_{in}}}}{{{R_{out}}}}} \right )\sqrt {{S_a}}. \end{equation}
Based on the mass balance,
$R_{in}$
in the droplet-radius-driven stage depends on its area due to the circle shape:
\begin{equation} {R_{in}} = \sqrt {\left ( {1 - {S_a}} \right ){A / \unicode{x03C0} }}. \end{equation}
Combining (3.6) and (3.13) in (3.12), we get the function of
$R^ *$
in the droplet-radius-driven stage:
\begin{equation} {R^ * } = \sqrt {{S_a}} \left ( {1 - \frac {{\sqrt {1 - {S_a}} }}{{\sqrt {{S_a}} }}} \right ) = \sqrt {{S_a}} - \sqrt {1 - {S_a}}. \end{equation}
By summing the release dynamics of trapped non-aqueous fluid to describe the evolution of
$R^*$
, we finally find that the release dynamics of trapped non-aqueous fluid is given by
\begin{equation} {R^ * } = \left \{ {\begin{array}{ll} {\sqrt {{S_a}} }&{( {0 \leqslant {S_a} \lt {A^ * }} )},\\ {\sqrt {{S_a}} - \sqrt {{\unicode{x03C0} / {{F_n}}} ( {{S_a} - {A^ * }} )} }&{( {{A^ * } \leqslant {S_a} \lt {A^ * } + {{{F_n}} / {2\unicode{x03C0} }}} )},\\ {\sqrt {{S_a}} - \sqrt {1 - {S_a}} }&{ ( {{A^ * } + {{{F_n}} / {2\unicode{x03C0} }} \leqslant {S_a} \leqslant 1} )}. \end{array}} \right . \end{equation}
The pore-shape-dependent angularity factor
$F_n$
and microscopic rough area factor
$A^*$
are the key geometrical parameters controlling the release process. Equation (3.15) can be used to describe the whole release process, including three stages: the film-driven, angular-radius-driven and droplet-radius-driven stages.
To validate these derived mathematical models, we conducted various microfluidic experiments and SC-LBM simulations with different pore geometries; detailed structure information is presented in Appendix B. During the experiments, all listed assumptions were observed to be valid, and all results are shown in Appendix B and supplementary movies 1–6. To quantify the self-removal processes in microfluidic experiments, the fluid–fluid interface evolution in the square (figure 7
a) and circle (figure 7
b) dead-end pores as typical dead-end structures were measured and compared with the mathematical model. All experimental data from dead-ends with varying pore shapes and microscopic rough areas collapse onto the theoretical prediction curves. Increasing microscopic rough area ratio
$A^*$
enhances self-removal outcomes (figures 7
a,b). Moreover, the release evolution curves differ significantly across dead ends with different pore shapes. For higher angularity factors
$F_n$
, detachment of the fluid–fluid interface from the microscopic structure occurs earlier, facilitating the self-removal process (figure 7
c). Figure 7(c) shows a more comprehensive prediction of the critical release evolution curves of different dead-end structures, where the self-removal process will go through three stages, classified as the film-driven, angular-radius-driven and droplet-radius-driven stages, which are described by (3.15). Consequently, the microscopic rough area ratio
$A^*$
and angularity
$F_n$
are the key parameters that control the capillary pressure gradient direction to trap or release the trapped fluid (figure 7
d).
Evolution curves and phase diagram of the transition from capillary trapping to complete release. (a,b) Prediction and fitting of different data sets, where the evolution curves of various experiments collapse onto the predicted curve in (a) the square and (b) the circle dead end (for the square,
$A^*= 0.21, 0.45, 0.58$
; for the circle,
$A^*= 0.32, 0.53, 0.60$
). The theoretical curves are obtained from (3.15). (c) The critical fluid release process for dead-end geometries with different shape factors and corresponding three stages of release processes. (d) Phase diagram of the transition from trapping to release from the dead end. The open and solid symbols represent the simulation (star-shaped symbols) and experimental cases (from triangle to circle symbols) where non-aqueous fluid was trapped or released from the dead end. The theoretical critical line is obtained from (3.16). Detailed geometry information about these dead ends is presented in figure 12, and the corresponding dynamic results of these dead ends are presented in figures 13–17 and the supplementary movies.
We summarize our experimental and simulation results using the phase diagram shown in figure 7(d). The open and solid symbols represent the experimental and simulation cases where non-aqueous fluid was trapped or released from the dead end. The circle symbols represent the microfluidic experimental results, and the star-shaped symbols represent the direct numerical simulation results. In all cases, we observe that the trapping state occurs at small values of
$A^*$
and
$F_n$
, as shown by the unfilled symbols in figure 3(e); however, as either
$A^*$
or
$F_n$
increases, complete release begins to be realized more easily, as shown by the solid symbols in figure 7(d). This is because increasing
$A^*$
will enhance the film-driven stage, and increasing
$F_n$
will enhance the angular-radius-driven stage. All the data appear to follow the predicted critical line given by (3.16), supporting our theory prediction, where
$A_{cri}^*$
is the critical microscopic rough area factor for releasing the trapped fluid from dead-end pores with different geometric shapes
$F_n$
:
\begin{equation} {A_{cri}^*} = \frac {{\unicode{x03C0} - {F_n}}}{{2\unicode{x03C0} }}. \end{equation}
This critical release line was derived from the continuous release evolution with
$R^*\geq0$
in (3.15). The excellent matching between the model predictions and experimental observations validates (3.16). According to the phase diagram in figure 7(d) and (3.16), there is an upper limit
${A_{cri}^*}={{{1 / 2} - {F_n}} / {2\unicode{x03C0} }}\lt 0.5$
. Based on this, some simple but useful conclusions can be made for engineering applications that when
${A^*} \gt 0.5$
, the non-aqueous phase can be released from the pores no matter what the shape of the dead end, while when
${A^*} \leqslant 0.5$
, it depends on the geometric shape, which can be determined based on figure 7(d) and (3.16).
3.4. Generality in complex porous structures
To determine the generality of the pore-scale results discussed thus far, we examined this anomalous capillary effect triggered by nanoparticle suspensions in random porous structures, in which dead-end geometries or dead-end-type porous structures are frequently and naturally occurring. Detailed information about the random porous structure has been shown in our previous work (Lei et al. Reference Lei, Liu, Xie, Yang, Wu and Wang2020). The porous media structure is regenerated based on statistical information about the sandstone obtained from the real geological system Changqing Oilfield in China. By colouring the fluids in the microfluidic porous media, we can distinguish invading and defending phases, and analyse the topology of each fluid. Therefore, the displacement efficiency and the topology of the defending fluid are described quantitatively by ganglion statistical analysis. The displacement experiments were performed on microchips with the same structure but different surface materials (silicon and silica) under a constant injection flow rate 1
$\unicode{x03BC} \text{l}\ \text{min}^-{^1}$
. Note that silica nanoparticles are easily adsorbed on the silica surface but not on the silicon surface.
Colour-modified image of the multiphase distribution and evolution of aqueous fluid saturation in (a,b) SiO
$_2$
porous media and (c,d) Si porous media. (e) Local view of nanoparticle-stimulated film flow and release phenomena in SiO
$_2$
porous media; positions 1 and 2 are chosen randomly from SiO
$_2$
-based microfluidic porous media.
The multiphase distributions at the final stage (1000 s) of nanoparticle suspension flooding and aqueous surfactant flooding in the SiO
$_2$
-based porous media are presented in figure 8(a), which indicates that nanoparticles can fragment the large oil cluster due to their enhanced film flow (figure 8
e). Figure 8(b) shows that the final displacement efficiency increases by approximately 19.6
$\%$
with nanoparticle suspension flooding compared to aqueous surfactant solution flooding. Similar experiments were conducted in Si-based porous media to elucidate the importance of nanoparticle adsorption, where silica nanoparticles cannot adsorb on the silicon surface (figure 8
c), and the final displacement efficiency increases by only approximately 7.5
$\%$
via the nanoparticle suspension (figure 8
d). This 7.5
$\%$
increase in displacement efficiency may be attributed to other factors related to the nanoparticles, such as particle adsorption kinetics at the fluid–fluid interface, and particle-ordering kinetics in the thin film (Lei et al. Reference Lei, Lu and Wang2023c
). However, this effect remains modest overall compared to the impact discussed in this paper. To further understand the pore-scale effect of non-aqueous fluid self-removal on increasing displacement efficiency, a local view picture was obtained from figure 8(a). The effects of the nanoparticle suspension and aqueous surfactant solution were compared: non-aqueous fluid removal via the nanoparticle suspension is the generic but essential phenomenon in porous media; nanoparticle-enhanced capillarity can promote breaking of large non-aqueous clusters into small areas, and fluid removal can frequently occur in dead-end or dead-end-type porous structures (figure 8
e). Figure 8(e) demonstrates directly that these results were within the expectation that frequent reversing capillary trapping in dead-end or dead-end-type porous structures originated from nanoparticle-stimulated film flow on the microstructure.
4. Conclusion
In summary, nanoparticle suspensions can trigger the release of non-aqueous fluid from dead-end pores. The multiscale capillarity on the hierarchical micro/nanostructures contributes directly to the reversing capillary trapping, as the nanostructure formed by nanoparticle adsorption can produce and maintain a robust hydrodynamic water film even on a microscopic convex surface under weakly water-wet conditions. With hydrodynamic film connectivity as a prerequisite, the capillary pressure gradient in the microstructure of the dead end drives the non-aqueous fluid motion out of the pore. The evolution dynamics of fluid release are described by (3.15), which is utterly theoretical, with no empirical parameters. The experimental data validate the developed mathematical models. Furthermore, a regime map is proposed to predict whether the trapped fluid can be entirely released, and all the data follow the critical release line (3.16). Finally, this nanoparticle-stimulated non-aqueous fluid release from a dead end is demonstrated as a generic phenomenon. Our findings offer insights into reversing capillary trapping during multiphase flow in the complex geometry, which can contribute to various engineering applications.
However, it should be noted that the formation of nanostructures on microscopic surfaces is closely related to the characteristics of nanoparticles, such as particle size, adsorption rate and concentration. In this paper, we did not analyse the dynamic process of nanoparticle adsorption for nanostructure formation, but only emphasized that the concave radius of the nanoparticle-formed nanostructure satisfies the capillary condensation condition
${r_{cc\_cri}} \lt {{2\sigma } / { ( {{P_c} - {\Pi ^{min }}} )}}$
. Based on the DLVO curve (figure 3
c) in this study, the concave curvature radius of the nanoscopic structure must be smaller than 200 nm; otherwise, this counterintuitive self-removal of non-aqueous fluid from the dead-end structure will not occur. The pore-scale model presented here is based on a polygon-shaped dead-end pore, so modifications to these models will be necessary when the dead-end pore shape is irregular. Under strong mainstream flow, snap-off phenomena frequently occur during these self-removal processes (Lei et al. Reference Lei, Lu, Liu and Wang2022, Reference Lei, Gong, Lu and Wang2023a
). The self-removal process driven by capillary force differences will need to be re-derived in such cases. These issues should be addressed in future work. Despite these considerations, our findings may provide valuable insights into multiphase flow manipulations for site-specific delivery, self-cleaning or self-removal. Furthermore, this research may offer a new perspective on solving scientific and technological challenges, such as soil and aquifer remediation, geological CO
$_2$
sequestration, hydrocarbon recovery, flow in fuel cell membranes, and hydration in biological systems, potentially leading to highly efficient fluid delivery and removal.
Supplementary movie.
Supplementary movie are available at https://doi.org/10.1017/jfm.2025.53.
Acknowledgements.
We extend our gratitude to Professor Wouter van der Wijngaart from KTH Royal Institute of Technology, Professor Chiyu Xie from Beihang University, Professor Wenbo Gong from China University of Petroleum, and Professor Jiangtao Zheng from China University of Mining and Technology (Beijing) for their valuable discussions and insights.
Funding.
We gratefully acknowledge support from the National Key Research and Development Program of China (no. 2019YFA0708704) and a grant from the National Science Foundation (NSF) of China (nos 12432013, 12272207).
Declaration of interests.
The authors report no conflict of interest.
Author contributions.
M.W. designed the research and supervised the project; W.L. performed research; X.L., G.Y. and S.B. analysed data; and W.L. and M.W. co-wrote the manuscript with feedback from all authors.
Appendix A. Validations of numerical solutions of the augmented Young–Laplace equation
The validity of the algorithm used to solve the augmented Young–Laplace equation was confirmed by analysing the thin-film distribution on pendular structures without considering extended DLVO forces (Sweeney et al. Reference Sweeney, Davis, Scriven and Zasadzinski1993), as well as the film thickness on spherical surfaces with extended DLVO forces. Theoretical results derived from the extended DLVO curves are presented in figure 3(c). Initially, we solved the Young–Laplace equation to estimate the thin-film configuration on pendular structures. In the case of triangular surface roughness, the contributions of pendular structures can be calculated analytically, and two solutions result:
\begin{equation} \left \{ {\begin{array}{*{20}{c}} {V_c^* = P_c^{* - 2}\left [ {m - \cos ^{ - 1}{{ ( {1 + {m^2}} )}^{ - {1 / 2}}}} \right ]},\\ {P_c^* \gt m{{ ( {1 + {m^2}} )}^{ - {1 / 2}}}}, \end{array}} \right . \end{equation}
\begin{equation} \left \{ {\begin{array}{*{20}{c}} {V_c^* = {{ ( {P_c^{* - 2} - 1} )}^{ - {1 / 2}}} + m - P_c^{* - 2}\sin^{ - 1}P_c^*},\\ {P_c^* \lt m{{ ( {1 + {m^2}} )}^{ - {1 / 2}}}}, \end{array}} \right . \end{equation}
where
$m$
is the ratio of the amplitude of surface roughness to half the wavelength of the roughness, dimensionless capillary pressure is defined as
$P_c^ * = {{{P_c}L} / \sigma }$
, dimensionless water film volume is defined as
$V_c^ * = {{{V_c}} / { ( {{L^2}w} )}}$
,
$V_c$
is the water film volume,
$L$
is a characteristic length scale of the solid geometry, and
$w$
is the surface width. The solutions at high values of capillary pressure correspond to menisci with a contact line lying along one of the flat surfaces. The contact line is free to move with changes in the capillary pressure. As the capillary pressure decreases, the contact line climbs up the surface until
$P_c^* = m{{ ( {1 + {m^2}} )}^{ - {1 / 2}}}$
. At this point, the contact line becomes pinned at the apex of the surface.
Validation of the algorithm for solving the augmented Young–Laplace equation. (a) The algorithm is validated against the analytical solution for periodic regular triangular surfaces provided by Sweeney et al. (Reference Sweeney, Davis, Scriven and Zasadzinski1993), without considering extended DLVO forces. (b) Thin-film thickness on spherical surfaces is analysed considering extended DLVO forces, where capillary pressure and spherical curvature radius are independent variables that control the film thickness.
Figure 9(a) compares the theoretical predictions with numerical results, where the lines represent theoretical predictions, and the points correspond to numerical results at various capillary pressure conditions. This demonstrates that when DLVO forces are negligible, the numerical algorithm for solving the Young–Laplace equation on complex structured solid surfaces is validated. By incorporating DLVO forces into the Young–Laplace equation, the thin-film thickness on spherical solid surfaces can be determined using capillary pressure and curvature radius through the augmented Young–Laplace equation. Figure 9(b) further compares the numerical solutions with theoretical results obtained from the extended DLVO curves, confirming the accuracy of the numerical solutions for the augmented Young–Laplace equation.
Appendix B. Multiphase lattice Boltzmann simulation
The SC-LBM was employed to simulate fluid–fluid interface evolution in a dead-end pore, which is widely used for modelling multiphase flow in microfluidics due to its high accuracy in capturing capillary phenomena (Shan & Chen Reference Shan and Chen1993; Zheng et al. Reference Zheng, Ju and Wang2018, Reference Zheng, Lei, Ju and Wang2021). For completeness, we briefly introduce the main framework of the model, with an emphasis on the forces of fluid and solid interaction. In the mesoscopic SC-LBM,
$k$
particle distribution functions are introduced to represent the multicomponent fluids. For this study, the multiple relaxation time collision model is utilized, which offers improved accuracy and stability by allowing the tuning of relaxation rates (Krüger et al. Reference Krüger, Kusumaatmaja, Kuzmin, Shardt, Silva and Viggen2017). The evolution function of
$f_i^k$
can be expressed as
\begin{equation} f_i^k ( {\boldsymbol {x} + \boldsymbol {e}_i\, \triangle {t} ,{t} + \triangle {t}} ) = f_i^k ( {\boldsymbol {x},t} ) - {\boldsymbol {M}^{ - 1}}\boldsymbol {S}\!\left (\! {\boldsymbol {M}\,{\boldsymbol {f}^k} ( {\boldsymbol {x},t} ) - {\boldsymbol {m}^{k,eq}}}\! \right ) + \triangle t\,{\boldsymbol {M}^{ - 1}}\!\left (\! {\boldsymbol {I} - \frac {\boldsymbol {S}}{2}} \!\right )\!{\hat {\boldsymbol {F}}^k}, \end{equation}
Two-phase simulation benchmark for the Young–Laplace equation. (a) Illustration of simulation set-up. The computational domain is cubic, and a bubble (red) with radius varying from 25 to 74 lattices is suspended in a fluid (blue). (b) The linear fitting curve verifies that the simulation follows the Young–Laplace equation. The slope corresponds to the dimensionless surface tension
$\sigma = 0.14$
.
Two-phase simulation benchmark for contact angle validation. (a) Schematic illustration and measurement of the contact angle. The sphere is fitted through interfacial points. (b) The linear fitting curve demonstrates that the contact angle remains consistent across different subcases, with the slope corresponding to contact angle
$\theta= 125^\circ$
. The computational domain is a cuboid, containing a droplet (red) with radius ranging from 20 to 61 lattices, attached to the wall within another fluid (blue).
Dead-end structures with different pore shapes (angularity factor
$F_n$
) and microscopic rough area (microscopic rough area ratio
$A^*$
). The critical release
$A_{cri}^*$
is determined by (3.16),
$A_{cri}^*= ({\unicode{x03C0} -F_n})/{2\unicode{x03C0} }$
.
where a diagonal matrix composed of relaxation rates of different moments
$\boldsymbol {S}$
can be expressed as
$\boldsymbol {S} = \text{diag} ( {S_\rho },{S_e},{S_{\varepsilon} },{S_j},{S_q},{S_j},{S_q},{S_j},{S_q},{S_\nu },{S_{\unicode{x03C0}} },{S_\nu },{S_{\unicode{x03C0}} },{S_\nu },{S_\nu },$
${S_\nu },{S_m},{S_m},{S_m} )$
. To achieve better stability and accuracy of multiphase flow and capillary phenomena, relaxation parameters are chosen as
$S_\rho=S_j=1$
,
$S_e= 1.19$
,
$S_{\unicode{x03C0}}=S_\varepsilon= 1.4$
,
$1/S_\nu={{3\nu }}/({{{c^2}\,\triangle t}})+1/2$
and
$S_m=S_q= {{8 ( {2 - {S_\nu }} )}}/({{8 - {S_\nu }})}$
(Chen et al. Reference Chen, Li, Valocchi and Christensen2017), where the lattice speed
$c = {{\triangle x} / {\triangle t}}$
is the ratio of lattice spacing
$\triangle x$
and time step
$\triangle t$
. Here,
$S_\nu$
is related to the kinematic viscosity
$\nu$
,
$\boldsymbol {M}$
is the transformation matrix that maps the distribution function to its moments,
${\boldsymbol {m}}^{k, eq}$
is a vector composed of equilibrium moments, and
${\hat {\boldsymbol {F}}^k} = \boldsymbol {M}{{\boldsymbol {F}}^k}$
accounts for body force in moment space (Guo, Zheng & Shi Reference Guo, Zheng and Shi2002). With the D3Q19 lattice,
$\boldsymbol {m}_i^{k,eq}$
and
$\boldsymbol {\hat {F}}_i^k$
are described as
\begin{equation} {\boldsymbol {m}^{k,eq}} = \left [\! {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {{\rho ^k}}\\ { - 11{\rho ^k} + 19{\rho ^k}( {u_x^2 + u_y^2 + u_z^2} )}\\ {\begin{array}{*{20}{c}} {3{\rho ^k} - \tfrac {{11}}{2}{\rho ^k}( {u_x^2 + u_y^2 + u_z^2} )}\\ {{\rho ^k}{u_x}}\\ {\begin{array}{*{20}{c}} { - \tfrac {2}{3}{\rho ^k}{u_x}}\\ {{\rho ^k}{u_y}}\\ {\begin{array}{*{20}{c}} { - \tfrac {2}{3}{\rho ^k}{u_y}}\\ {{\rho ^k}{u_z}}\\ {\begin{array}{*{20}{c}} { - \tfrac {2}{3}{\rho ^k}{u_z}}\\ {{\rho ^k}( {2u_x^2 - u_y^2 - u_z^2} )}\\ { - \tfrac {1}{2}{\rho ^k}( {2u_x^2 - u_y^2 - u_z^2} )} \end{array}} \end{array}} \end{array}} \end{array}} \end{array}}\\ {{\rho ^k}( {u_y^2 - u_z^2} )}\\ {\begin{array}{*{20}{c}} { - \tfrac {1}{2}{\rho ^k}( {u_y^2 - u_z^2} )}\\ {{\rho ^k}{u_x}{u_y}}\\ {\begin{array}{*{20}{c}} {{\rho ^k}{u_y}{u_z}}\\ {{\rho ^k}{u_x}{u_z}}\\ {\begin{array}{*{20}{c}} 0\\ 0\\ 0 \end{array}} \end{array}} \end{array}} \end{array}}\! \right ],\; {{\hat {\boldsymbol {F}}}^k} = \left [\! {\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} 0\\ {38\boldsymbol {u} \boldsymbol{\cdot} \boldsymbol {F}}\\ {\begin{array}{*{20}{c}} { - 11\boldsymbol {u} \boldsymbol{\cdot} \boldsymbol {F}}\\ {{F_x}}\\ {\begin{array}{*{20}{c}} { - \tfrac {2}{3}{F_x}}\\ {{F_y}}\\ {\begin{array}{*{20}{c}} { - \tfrac {2}{3}{F_y}}\\ {{F_z}}\\ {\begin{array}{*{20}{c}} { - \tfrac {2}{3}{F_z}}\\ {2( {2{F_x}{u_x} - {F_y}{u_y} - {F_z}{u_z}} )}\\ { - ( {2{F_x}{u_x} - {F_y}{u_y} - {F_z}{u_z}} )} \end{array}} \end{array}} \end{array}} \end{array}} \end{array}}\\ {2( {{F_y}{u_y} - {F_z}{u_z}} )}\\ {\begin{array}{*{20}{c}} { - ( {{F_y}{u_y} - {F_z}{u_z}} )}\\ {{F_x}{u_x} + {F_y}{u_y}}\\ {\begin{array}{*{20}{c}} {{F_y}{u_y} + {F_z}{u_z}}\\ {{F_x}{u_x} + {F_z}{u_z}}\\ {\begin{array}{*{20}{c}} 0\\ 0\\ 0 \end{array}} \end{array}} \end{array}} \end{array}}\! \right ]\!, \end{equation}
where
$\boldsymbol {F}=[F_x, F_y, F_z]$
is the body forcing term,
$\boldsymbol {u}=[u_x, u_y, u_z]$
is the velocity, and
$\rho$
is the density.
The release process from the triangle dead-end structure via nanoparticle suspension, compared with the aqueous surfactant solution effect.
The release process from the square dead-end structure via nanoparticle suspension, compared with the aqueous surfactant solution effect.
The release process from the pentagon dead-end structure via nanoparticle suspension, compared with the aqueous surfactant solution effect.
The release process from the octagon dead-end structure via nanoparticle suspension, compared with the aqueous surfactant solution effect.
The release process from the circle dead-end structure via nanoparticle suspension, compared with the aqueous surfactant solution effect.
An additional fluid–fluid interaction force was introduced apart from other body forces to achieve phase separation in SC-LBM. Since there are only two components in our system, the fluid–fluid force can be expressed as
\begin{equation} \boldsymbol {F}^k = - {\psi ^k} ( {\boldsymbol x} )\sum _i {{G_i}\,{\psi ^{\bar k}} ( {\boldsymbol x + {{\boldsymbol {e}_i}}\, \triangle {t}} )} \,{{\boldsymbol {e}_i}}, \end{equation}
where
$\psi$
is the pseudo-potential representing the effective number density of the component
$k$
, and the other component
$\overline k$
. In this work,
$\psi =\rho$
is chosen for simplicity, and
$G_i$
is the coefficient that relates to the surface tension,
\begin{equation} {G_i} ( {{\boldsymbol x},{{\boldsymbol x}^\prime }} ) = \left \{ {\begin{array}{ll} {{w_i}G}, &{{{\boldsymbol x}^\prime } = {\boldsymbol x} + {{\boldsymbol e}_i} },\\ 0, &\text{otherwise}, \end{array}} \right . \end{equation}
where
$w_i$
is the lattice weight. As a result, the fluid–fluid interaction force is a short-range force that only affects the interfacial area due to the density difference.
The related macroscopic variables for density
$\rho$
, momentum
$\rho \boldsymbol {u}$
and total velocity
$\boldsymbol {u}$
are defined as
\begin{equation} {\rho ^k} = \sum {f_i^k} ,\quad{\rho ^k}{\boldsymbol {u}^k} = \sum _i {f_i^k {{\boldsymbol {e}_i}}} + \frac {{{\rho ^k}{\boldsymbol {F}^k}}\,\triangle t}{2}, \end{equation}
\begin{equation} \rho = \sum \limits _k {{\rho ^k}} ,\quad\rho \boldsymbol {u} = \sum \limits _k {{\rho ^k}{\boldsymbol {u}^k}}. \end{equation}
Additionally, the virtual fluid method is employed to realize the wettability of the system. This method assigns virtual fluid density at solid boundaries according to the average of the surrounding
$k$
th component density:
\begin{equation} \left \{ {\begin{array}{*{20}{c}} {\rho _{avg}^k = \dfrac {1}{{NP}}\sum \rho^k \left ( {\boldsymbol {x} + {{\boldsymbol {e}}_{{i}}}\, \triangle {t}} \right )},\\[11pt] {\rho _s^k = \rho _{avg}^k + \triangle {S^k}},\\[7pt] {\triangle {S^w} = - \triangle {S^{nw}} = \triangle S}, \end{array}} \right . \end{equation}
where
$NP$
is the number of fluid nodes that directly connect to the solid node,
$nw$
denotes the non-wetting phase, and
$w$
denotes the wetting phase. The density of the non-wetting phase becomes smaller, while the density of the wetting phase becomes larger at the solid boundary node. As a result, the solid surface is more affinitive with the wetting phase.
Several benchmarks were conducted to validate our code, including Young–Laplace and contact angle validations. The first benchmark involved calculating the relationship between droplet radius and pressure difference to validate the Young–Laplace equation. The physical model for this benchmark is shown in figure 10. In this model, a spherical droplet (red) is placed at the centre of another fluid (blue). At equilibrium, the pressure difference inside and outside the droplet, which can also be obtained through simulation, follows the Young–Laplace equation:
\begin{equation} \triangle p=\sigma \left (\frac {1}{R_1}+\frac {1}{R_2}\right )=\frac {2\sigma }{R}, \end{equation}
where
$\sigma$
is surface tension, and
$R$
is the radius of the bubble. A cubic domain of
$240\times240\times240$
lattices with a periodic boundary condition at all cubic boundaries is utilized in this simulation. According to previous analyses, the significant phenomenon of this simulation is a quasi-static process driven by capillary force. Thus the absolute values of density and viscosity are insignificant. In this simulation, the density of the two phases is 1.0, and kinetic viscosity is
$\nu=0.1667$
(normalized by
${\triangle {x^2}} / {\triangle t}$
).
$G = 5.0$
is chosen to separate the two phases. Five subcases with radii 25, 37, 49, 62 and 74 lattices are calculated correspondingly.
To verify the accuracy of the wettability representation used in the simulation, we conduct the second benchmark for contact angle validation. The schematic of this benchmark calculation is shown in figure 11, where a droplet with a different radius is attached to a flat wall. As a virtual density method, a predetermined
$\triangle S$
corresponds to a certain contact angle, which models a certain two-phase fluid–solid system. Thus the contact angle
$\theta$
can be computed directly from simulations at a steady state by the geometrical relations
\begin{equation} \cos \theta = - \frac {{{y_c}}}{R}. \end{equation}
Since gravity is neglected in the benchmark, the droplet assumes a spherical shape. By extracting all interfacial points, the radius and centre of the sphere can be determined using least squares fitting. The
$y$
-coordinate of the sphere’s centre is denoted as
$y_c$
, and the radius is denoted as
$R$
. A cuboid domain of
$240\times240\times240$
lattices is used. The viscosity and density are identical to those in the Young–Laplace benchmark. A value
$\triangle S = 0.1$
is selected to achieve the desired contact angle.
Lattice Boltzmann simulations about the release process of non-aqueous fluid from the triangle and square dead-end structures are presented in movie 6 of the supplementary material. The self-removal results of these two simulation cases have been presented in the phase diagram as star-shaped symbols in figure 7(d).
Appendix C. Dead-end structure design and detailed experimental results of releasing non-aqueous fluid from various dead-end structures
Various dead-end microchips were designed to observe the self-removal process, as shown in figure 12. The results based on these dead-end structures define the parameter space
${F_n}$
–
$A^*$
, as presented in figure 7(d). Microfluidic experiments demonstrate the self-removal process of non-aqueous fluids from these structures, which are detailed in the supplementary materials (triangle in movie 1, square in movie 2, pentagon in movie 3, octagon in movie 4, and circle in movie 5), using (a) nanoparticle suspensions and (b) aqueous surfactant solutions. Representative images of the fluid–fluid interface at various stages of the experiment are summarized in figures 13–17. For all conditions (
$F_n$
and
$A^*$
), aqueous surfactant solutions had no effect, while nanoparticle suspensions promoted film spreading within the dead-end structures. The final self-removal outcomes are governed by
$F_n$
and
$A^*$
. All self-removal processes follow three typical stages, as introduced in § 3.3. However, because the circular dead-end lacks angularity, the process skips the angular-radius-driven stage (figure 17). As the dead-end structure transitions from triangle to circle, the critical microscopic rough area ratio
$A_{cri}^*$
increases, aligning well with our theoretical predictions.
