A many-body dissipative particle dynamics study of fluid–fluid spontaneous capillary displacement

Abstract

The fluid–fluid spontaneous capillary displacement is a model related to many important natural processes and industrial applications, such as groundwater remediation and oil recovery. To study the displacement process in microscale where experimental observations are hard to perform, computational simulations are an adequate approach. In the present work, we report a many-body dissipative particle dynamics (MDPD) simulation study on the fluid–fluid spontaneous capillary displacement. System parameters such as the wetting property of both fluids, the miscibility between them, and the capillary radius, can be adjusted independently and separately in the simulation, thus their influences on the spontaneous capillary displacement process can be observed individually. Based on our previous study of MDPD simulations on fluid imbibition, we are able to monitor the motion of the fluid–fluid interface in this spontaneous capillary displacement simulation accurately. Combined with the new method for three-phase contact angle calculation proposed in the present study, the displacement curve obtained is in good accordance with the theoretical description even in the very beginning of the spontaneous capillary displacement motion. Our findings provide microscopic information about the fluid–fluid spontaneous displacement process, and lay a basis for further studies on more complex systems.

---

Introduction

Fluid–fluid displacement is a process taking place in a capillary that involves simultaneously the imbibition of the displacing phase (denoted as fluid A or FA) and the drainage of the displaced phase (denoted as fluid B or FB). FA flows into the capillary continuously from a reservoir, while FB flows out at the same time into another, which leads to an overall effect that FB in the capillary is replaced by FA.

The fluid–fluid displacement can roughly be classified into spontaneous displacement or forced displacement, depending on whether an external force is applied to trigger the process. From the standpoint of energy, if the capillary inner surface covered with FA is of lower system energy than that covered with FB, the pre-occupying FB will be replaced by FA, i.e., the displacement process would take place spontaneously. On the contrary, if the inner surface of capillary favors FB, extra energy or external force would be required to initiate the displacement.

The fluid–fluid displacement in capillary is closely related to many natural and industrial processes. For instance, in groundwater remediation, the ultimate target is a complete removal of the non-aqueous phase liquids (NAPLs), a class of hazardous organics immiscible with water and sticking in the porous aquifer.^1,2 Similarly, in oil recovery, petroleum deeply trapped in the fractured field is expected to be entirely extracted.^3,4 The common feature of both processes is that the non-aqueous phase is replaced by water or aqueous solution. These processes have been attracting considerable attentions of the scientific and engineering community,^5–12 because, for such important environment or energy issues, any improvement in the efficiency or yield is of utmost significance, and the fundamental understanding of the displacement mechanism, in particular in microscale, is of high value.

The behavior of fluid–fluid displacement is affected by various properties, including the wettability of the FA–capillary and FB–capillary interfaces, the miscibility of the FA–FB interface, the viscosity of both fluids, and the radius of the capillary, etc.^13,14 In most cases of practical experiments, independently varying these properties would not be an easy task. For example, highly efficient surfactants are usually dissolved into the solution to manipulate the fluid–fluid interfacial tension.^15–19 However the surfactant could also be very likely to influence the FA–capillary interface, and to alter the wetting properties of the displacing solution. Generally speaking, these influential elements are entangled with each other, and the individual characters remain obscure in experimental studies.

Moreover, although it is possible today to perform observations on experimental research in fine channels at micrometer scale, precise measurements for displacement processes in even thinner slit (e.g., tens or hundreds of nanometers) remains highly challenging.^20,21 Information for flow patterns at this scale are rather crucial for applications such as oil recovery in ultra-low permeability oil reservoirs, which account for a significant sector of oil exploration in China. In view of this, the computational simulation is an adequate approach,^22–26 in which any independent variation in related properties would be quite natural and straightforward, and it permits the fluid–fluid displacement process in capillary to be observed in molecular details.

Previously, we have conducted simulations on forced fluid–fluid displacement in capillary.²⁷ The roles played by several influential elements were analyzed, and their effects on the starting force and the residual oil content have been discussed in detail. However, in most cases of the forced displacement, due to the presence of relatively large external pressure, the interfaces therein deviates a lot from a spherical shape,²⁸ which makes the simulated displacement process very difficult to be compared with the theoretical description using known formulae. In the case of spontaneous capillary displacement, the absent of external pressure permits the observation of the evolution and dynamics of the meniscus, which thus brings the possibility for deeper analysis.

In the present paper, we report a systematic simulation study on spontaneous capillary displacement, where the capillary employed have a round-shaped slit. Several fluid properties in the process are varied independently to collect information about the displacement under each condition. The spontaneous capillary displacement process is further analyzed theoretically based on an expression derived by extending single-phase fluid motions in capillary^29–32 to the two-fluid one, and it turns out that the simulations are well consistent with the theoretical descriptions.

Theory

Before analyzing the process of spontaneous capillary displacement, two basic fluid motions in capillary should be mentioned first, i.e., the spontaneous capillary imbibition (SCI) and the spontaneous capillary drainage (SCD). As we know, it is the wetting property that determines the dynamic process of a specific fluid.³³ Generally speaking, when a fluid forms a static contact angle (θ₀) less than 90°, the fluid is usually said to be wetting for the capillary, otherwise the fluid is regarded as non-wetting.³⁴ Note that the term “completely wetting” corresponds to the case of θ₀ = 0°, while other wetting cases are referred as “partially wetting”.

The SCI process takes place in the capillary only for wetting fluids, while the SCD happens only for non-wetting fluids. When an empty capillary is brought in touch with a reservoir containing a wetting fluid, the fluid will flow spontaneous into the capillary, which is referred as SCI. On the contrary, when a capillary filled with non-wetting fluid is brought in touch with its reservoir, the fluid will flow continuously out of the capillary, which is referred as SCD. In both cases, the driving capillary force (ΔP) can readily be calculated based on the radius of capillary (r), the surface tension of the fluid (γ), and the contact angle (θ) through the Young–Laplace equation:³³

The resistance in both SCI and SCD has two sources, i.e., the viscous one and the inertial one. The viscous resistance stems from the friction between contacting fluid layers in capillary, which can be calculated according to the Hagen–Poiseuille equation,³⁵ and for completely wetting case P_Viscous is:

where η is the viscosity, l_c the fluid length in capillary. The inertial resistance is generated due to the acceleration of the fluid and can be calculated according to the following formula:^29,30where dp(t)/dt is the instantaneous momentum change, which consists of two parts:³⁶where l_c and l_r are fluid length in capillary and in reservoir, respectively, while m_c and m_r are the corresponding fluid mass. The viscous and inertial resistances contribute differently during different stages of the SCI and SCD processes. The viscous resistance increases or decreases steadily with increasing or decreasing the length of fluid column in SCI or SCD,³⁰ while in contrast, the inertial resistance only contributes significantly at the very beginning stage of the SCI and SCD processes owing to the dramatic velocity change at that moment. For experiments and long-time simulations, the inertial resistance could be neglected, but recent studies have revealed that it is indispensable for precise description of the fluid motion at the initial stage.³⁶

A balance between the driving force and the resistance leads to a generalized formula for SCI and SCD (in differential equation form):

By further assuming the inertial contribution is negligible, eqn (5) changes into:

When solving analytically for meniscus height l_c as a function of time t, one can arrive at the famous Lucas–Washburn equation^31,32 (eqn (7)), which states that l_c is proportional to the square root of t.

The Lucas–Washburn equation has been widely employed in describing the SCI process for completely wetting fluids. The validity of this equation has also been verified in various simulation studies.^{29,30,36–41} However, Supple and Quirke found that the height of oil imbibed in a carbon nanotube was not proportional to the square root of time, but obeyed a linear relation with time.^42,43 They have attributed this result to the short time scale and the ultra-low friction of graphitic surface, and thought that for rougher surfaces a Washburn-like equation would be obeyed.

Note that, in simulation studies, several researchers have suggested and demonstrated that, the contact angle θ in Lucas–Washburn equation should be the dynamic contact angle (θ_d) rather than θ₀ that usually employed in experiments, so that the instantaneous variation in driving force can be rationally included.³⁷ Besides, one of the boundary conditions chosen to solve the differential equation is so-called “no-slip” condition, which states that the fluid will have zero velocity at the fluid–capillary interface. This condition is correct only for completely wetting fluids but not for other wetting cases. To account for the fluid slip, a term “slip length” b has been introduced, and the thus-derived equation for both SCI and SCD should reads:^30,32

We have shown in previous simulation study that,³⁰ with appropriate slip length chosen, the fluid motions in SCI and SCD can be correctly described by the above equation.

The success of SCI and SCD simulations brings possibility for a systematic simulation study of the spontaneous capillary displacement. In certain sense, the spontaneous capillary displacement can be regarded as a combination of SCI and SCD. However, for spontaneous capillary displacement, when FA and FB coexist at equilibrium in the capillary and are closely attached to each other, there is no longer static contact angle for either FA or FB, but rather a three-phase static contact angle (θ₀₃) formed at the FA–FB–capillary three-phase interface. Note that, throughout this paper, θ₀₃ is defined according to the displacing phase (FA), thus if θ₀₃ is less than 90°, FA will displace FB spontaneously. Therefore, the capillary pressure should be calculated with the fluid–fluid interfacial tension γ_AB and the three-phase dynamic contact angle (θ_d3) using the Young–Laplace equation (eqn (1)). As that in the SCI and SCD processes, the resistance in spontaneous capillary displacement stems from the viscous and the inertial resistances, which are contributed by both FA and FB. In the present work, both FA and FB are set to have the same viscosity for simplicity. The case for FA and FB with different viscosities is left for later investigation. Considering θ_d3 and the slip of fluids, the differential equation for fluid–fluid spontaneous capillary displacement writes:

where b_A and b_B are slip lengths for FA and FB, and dp(t)/dt is calculated following:where c and r in the subscript refers to capillary and reservoir, while A and B refers to contributions from FA and FB, respectively.

If similar assumptions are applied as that used in deriving Lucas–Washburn equation, i.e. neglecting the inertial contribution and considering the θ_d3 kept almost unchanged from a constant value θ, eqn (9) then turns into the following equation:

In this formula, the moving speeds for FA and FB are equal, i.e. dl_cA/dt = −dl_cB/dt. If we further assume that (r² + 4b_Ar)⁻¹ ≈ (r² + 4b_Br)⁻¹ = (r² + 4br)⁻¹, then it could be reduced to:

Since the sum of l_cA and l_cB remains constant as the length of the capillary l_c, integrating this formula one gets

This expression shows a linear dependence of l_cA on t. As we will show in the later part of the present work, under steady state of the spontaneous capillary displacement process, the θ_d3 will fluctuate around such constant value, which validates our assumption here. Moreover, this linear dependence of l_cA on t is also evident from the recorded displacing curves along simulation beyond the very beginning accelerating stage.

During the simulation, l_cA is recorded every 0.1 MDPD time unit by monitoring the volume equivalent front of FA–FB interface,³⁶ and according to the constant amount of FA and FB particles as well as known density, other fluid length and mass either in capillary or reservoirs are readily calculated, which permits the momentum changed at each moment to be determined. Additionally, the θ_d3 is extracted along the process and later substituted into eqn (9) at each time point. eqnn (9) is a second-order ordinary differential equation taking l_cA and t as the dependent and independent variables, respectively. The equation could be solved numerically using computer for t spans 0 to 400 with two initial conditions: (1) l_cW = l_0cA at t = 0, where l_0cA is the initial length of water column in capillary before displacement starts; (2) dl_cA/dt = 0 at t = 0, which means the initial velocity is zero. The detailed comparison of simulation results and the numerical calculation of eqn (9) will be given in the Results and discussion section.

Computational method

Many-body dissipative particle dynamics (MDPD)

As originally proposed by Hoogerbrugge and Koelman,⁴⁴ and later advanced by Groot and Warren⁴⁵ in elucidating the underlying physical essence, the dissipative particle dynamics (DPD) method is versatile in molecular simulations. Each DPD particles is characterized by its position r_i, velocity v_i, and mass m_i, and the movement is governed by Newton's second law. The forces that act on a DPD particle can be divided into conservative, dissipative and random parts, i.e., F^C_ij, F^D_ij, and F^R_ij, respectively, with a cutoff distance r_c (reduced unit of length) imposed on all force components:

F_ij = F^C_ij + F^D_ij + F^R_ij (14)

The conservative part is designed as a soft repulsive force,

F^C_ij = A_ijw_C(r_ij)e_ij (15)

where A_ij is the maximum repulsion between particle i and j, w_C(r_ij) is a weight function vanishing at r > r_c = 1, and r_ij = r_j − r_i, r_ij = |r_ij|, e_ij = r_ij/r_ij. The dissipative part and random part write as follows:

F^D_ij = −γw_D(r_ij)(v_ij·e_ij)e_ij (16)

F^R_ij = σw_R(r_ij)ξ_ije_ij (17)

where the coefficient γ describes the extent of dissipation, ξ_ij is a random number obeying Gaussian distribution with the mean value, amplitude and variance being 0, σ and 1/δt, respectively, w_D(r_ij) and w_R(r_ij) are weight functions vanishing beyond cutoff distance r_c. It had been shown by Español and Warren⁴⁶ that the dissipative and random force related to each other through w_D(r_ij) = [w_C(r_ij)]² and σ² = 2γk_BT, such that a thermostat is constituted that conserves the momentum and does not influence the system equilibrium thermodynamics.

Although the DPD method is efficient and has been widely applied to many practical problems, Louis et al.⁴⁷ had pointed out that the quadratic equation of state (EOS) derived from the conservative force does not contain a van der Waals loop, which fails to represent the EOS of real fluid. The dissipative and random forces together act as a thermostat, which leaves the effective force between DPD particles solely repulsive, leading to the simulation box filled fully with particles due to the lacking of force making them stick together. Consequently, the original DPD method is not qualified for studying systems with free surfaces (e.g. liquid–gas interface).

Generally speaking, any attempt that tends to equip DPD with a capability by incorporating a many-body force is called MDPD. Warren's MDPD extension⁴⁸ has been well justified and easy to be implemented, thus we employed it in this work. In his model, while preserving the F^D_ij and F^R_ij force components, the conservative part is replaced with:

F^C_ij = A_ijw_C(r_ij)e_ij + B([small rho, Greek, macron]_i + [small rho, Greek, macron]_j)w_d(r_ij)e_ij (18)

where w_C(r_ij) and w_d(r_ij) are weight functions vanishing at distance beyond cutoff r_C and r_d, respectively.

w_C(r) = 1 − r/r_c (19)

w_d(r) = 1 − r/r_d (20)

The partial amplitude ρ depends on a weighted local density and for the ith particle its value can be obtained as follows:

w_ρ(r) = 15(1 − r/r_d)²/(2πr_d³) (22)

Warren's MDPD method simply changes the sign of the original DPD's conservative force (the first term in eqn (18) with A_ij < 0, cutoff r_c = 1) and adds an additional short-ranged many-body repulsive force (the second term in eqn (18) with B > 0 and r_d = 0.75). In the present study, the timestep is set as 0.01 MDPD time unit, B is kept at 25 while A is set at −40 for interactions between the same kinds of particles, which lead to an equilibrium density of 6.00 ± 0.05. The surface tension and viscosity are determined to be γ = 7.62 ± 0.05 and η = 7.39 ± 0.06 in individual simulations following the literatures.^49,50 The tunable interaction parameter between solid–liquid are denoted as A_SL, which varies between −40 and 0. The interfacial tensions between two fluids are treated as that in our previous report for the forced capillary displacement.²⁷ All simulations were carried out using the LAMMPS package⁵¹ with MDPD implemented by modifying the DPD codes therein, and snapshots were taken in VMD.⁵²

Results and discussion

Interfacial tension

The fluid–fluid interfacial tension γ_AB is a key parameter determining the driving force for the fluid–fluid spontaneous capillary displacement according to the Young–Laplace equation (eqn (1)). Therefore, the calculation of interfacial tension is the first requirement for the correct description of this process. In our previous study on forced capillary displacement, we have demonstrated a method for the calculation of interfacial tensions.²⁷ Briefly, the simulation system consists of two fluid slabs (denote as FA & FB) sticking to each other and enough vacuum regions were left to each side of the fluid slabs. The calculated total tension γ_Total comprises three parts, i.e., the FA surface tension γ_A, the FB surface tension γ_B, and the FA–FB interfacial tension γ_AB, which are contributed from the leftmost, the right most, and the intermediate interfaces, respectively. Since the surface tensions of FA and FB are obtained from independent surface tension calculations of each single component (γ_A = γ_B = 7.62), the γ_AB can be obtained by subtracting γ_A and γ_B from the calculated γ_Total.

The thus-obtained γ_AB increases regularly (from 0 to ca. 15) with the fluid–fluid interaction parameter A_AB (from −40 to 0, representing a decline of attracting force between the FA and FB particles). The two limit choices of A_AB (−40 and 0) correspond to a fully miscible and an entirely immiscible configuration at the interface respectively. The FA and FB slabs would detach from each other when A_AB = 0, since there is not enough force to hold them together. Table S1 in the ESI† lists specific γ_AB values for several typical A_AB choices, which will later be substituted into eqn (9) as a constant.

Slip length

The second requirement for the correct description of the spontaneous capillary displacement process is to determine the slip length of fluid with different wettability, which has been achieved in our previous study of spontaneous capillary imbibition and drainage.³⁰ Although most published simulation research are focused on the motion of completely wetting fluid which can be well described with classic Lucas–Washburn equation, it was found that, only by including the slip length, can the thus-obtained eqn (8) be successfully used to describe the SCI or SCD process for completely/partially wetting or non-wetting fluids. This slip length is fitted during individual MDPD simulation for fluids with a specific wettability, it is small for completely wetting and partially wetting fluids, and increases dramatically for partially non-wetting fluids with increasing A_SL. Table S2 in the ESI† lists specific slip length values of fluid with typical interaction parameter A_SL, which will later be used as known constants to solve the differential equations for spontaneous capillary displacement process.

Effective fluid front

To appropriately determine the fluid front (i.e., the representing position of the interface) is the third requirement for the correct description of the spontaneous capillary displacement process. In our recent report regarding the effective fluid front of a moving meniscus in capillary, we have proposed a method to properly determine the representing position of the meniscus front.³⁶ In simulations, a fluid cylinder can be divided into concentric layers. We found that neither the outmost layer nor the central one gave a correct representation of the meniscus movement in the simulation of SCI, while a layer somewhere in between showed the best consistency with the theoretical description. The position of such a layer has been determined, in view of volume equivalent, to be at 72.1% of the radius away from the axis. This method has been extensively verified with MDPD simulations for SCI at different capillary radius.³⁶ For all simulations performed in this paper, the motion of the fluid–fluid interface in capillary is described by monitoring the effective fluid front of FA, i.e., the front point of a FA layer located at 72.1% of the radius away from the axis. The displacement curve is defined as the change of this effective fluid front as a function of time.

Three-phase static contact angle (θ₀₃)

The simulation of the spontaneous capillary displacement process needs the fluids to start from their static equilibrated state. This means that obtaining the three-phase static contact angle (θ₀₃) beforehand is another requirement for the spontaneous capillary displacement simulation. The θ₀₃ is formed at the three-phase interface (i.e., FA, FB, and capillary wall), whose magnitude is affected by the wetting properties of FA and FB and the miscibility between them.

The θ₀₃ can be obtained in the experiment by simply immersing the planar substrate and the attached droplet into another fluid bath,¹² based on which a simulation version could be constructed through imitation: in a simulation box, a flat substrate is partially covered with FB, and a small FA droplet is immersed into the FB and onto the substrate. After sufficient equilibration, the θ₀₃ could be extracted from this simulation. Albeit straightforward, this method has several shortcomings. First, it is computationally intensive, because a large simulation substrate is needed to avoid interactions across the periodic boundaries between the droplet and its image, and the most parts of the computational power is wasted at the large FB reservoir. Second, the significant mobility of FA droplet (especially in the partially non-wetting case) requires additional techniques to hold the droplet quasi-stationary, so that averaging between configurations is possible. Third, the θ₀₃ thus obtained may not be applicable to a capillary system, whose inner wall is curved.

Here we propose a novel method for the θ₀₃ calculation based on a capillary system. Since this part is not the main theme of this paper, we give a brief description here. The initial state of the simulation box consists of five components (Fig. 1): (1) a capillary with cylinder inner shape and rectangular outer wall; (2) a wall placed at the left entrance of the capillary; (3) a FA cylinder filled in the left side of the capillary and attached to the wall; (4) a FB cylinder filled in the capillary and attached to the right side of the FA cylinder; (5) a vertical and rigid piston attached to the right side of the FB cylinder and free to move along the axis in the capillary. Note that the freely moving, vertical piston in this configuration is necessary, such that the right interface of the FB cylinder will remain vertical during the simulation and will not impact on the θ₀₃ calculation.

Fig. 1 Simulation system for the calculation of three-phase static contact angle (θ₀₃). (a) Pre-equilibrated homogeneous system ready for partition. (b) Cross-section view of the capillary. (c) and (d) are the slide view of the system at the initial and final state of simulation, respectively. Particles are divided into five types: a capillary (black) with a round slit; a wall (green) sealing the left entrance of the capillary; a FA column (iceblue) attached to the wall; an FB column (brown) attached to the FA; a free piston (purple) attached to the right-hand side of the FB.

Table 1 lists the thus-obtained θ₀₃ values for different types of FA and FB. In these calculations, while the FA–wall and FB–piston interaction parameter is set as constants A_AW = A_BP = −40, the FA–capillary interaction A_A, FB–capillary interaction A_B, and FA–FB interaction A_AB are adjusted to mimic different properties of the displacing and displaced phase. It is notable that the fluid–fluid interaction A_AB plays an important role in determining the θ₀₃. In extreme cases of A_AB = 0 or −40, the FA–FB interface will break into two separate surfaces or completely disappear, thus the so-called interface no longer exists. Therefore the θ₀₃ values here listed are calculated for interfaces with A_AB being −20.

Table 1 Three-phase static contact angle (θ₀₃) calculated with different sets of parameters

Series 1^a		Series 2^b		Series 3^c		Series 4^d
r	θ03 (°)	AA	θ03 (°)	AB	θ03 (°)	AAB	θ03 (°)
a For calculations in this series, parameters are set as: AA = −40, AB = −30, AAB = −20.b For calculations in this series, parameters are set as: AB = −25, AAB = −20, r = 10.c For calculations in this series, parameters are set as: AA = −40, AAB = −20, r = 10.d For calculations in this series, parameters are set as: AA = −40, AB = −30, r = 10.
5	39.00 ± 0.93	−40	20.06 ± 1.03	−37.5	76.32 ± 1.25	−27.5	14.99 ± 1.89
10	42.78 ± 1.42	−37.5	41.29 ± 0.89	−35	64.00 ± 0.69	−25	28.79 ± 1.61
15	44.32 ± 0.71	−35	55.61 ± 1.91	−32.5	53.21 ± 1.26	−22.5	36.09 ± 1.50
20	43.40 ± 1.03	−32.5	65.70 ± 1.14	−30	42.78 ± 1.42	−20	42.78 ± 1.42
25	43.98 ± 0.45	−30	75.81 ± 1.32	−27.5	30.21 ± 0.26	−17.5	45.47 ± 1.10
		−27.5	83.28 ± 1.07	−25	20.06 ± 1.03	−15	47.33 ± 0.85

---

Spontaneous capillary displacement

The simulation model for spontaneous capillary displacement consists of 3 components, as shown in Fig. 2: (1) a capillary of square cross-section with a round-shaped slit; (2) a FA reservoir connected to the left entrance of the capillary; (3) an FB reservoir connected to the right exit of the capillary, with the capillary being pre-filled with FB. The system construction process is detailed in the ESI.†

Fig. 2 MDPD simulations for the spontaneous capillary displacement. (a) The system consists of a capillary (in 40 MDPD length unit) prefilled with FB, a FA reservoir (in 30 MDPD length unit), and an FB reservoir (in 10 MDPD length unit). A θ₀₃ is formed with the arc bottom aligned to the capillary left entrance before the spontaneous capillary displacement starts. Capillary, FA, and FB particles are colored black, iceblue, and brown, respectively. (b), (c) and (d) are system states at t = 0, t = 200, and t = 400, respectively. (e) Cross-section view of the vacuum capillary. Displacing curve is recorded by monitoring the effective fluid front point (i.e., a position at the FA–FB interface, whose radial coordinate is 0.721),³⁶ which is marked as red points in (a) to (d).

As previously discussed, to trigger the spontaneous capillary displacement process, FA should have stronger tendency than FB to spread on the inner surface of the capillary. To describe this, we let A_A be always less than A_B. A typical spontaneous capillary displacement process can be divided into two stages: Stage 1 ranges from t = 0 to t = ∼80 (for capillary of r = 10), and is the acceleration process of fluids; Stage 2 covers the time behind, when the fluids flow at a quasi-constant speed. These two stages can be clearly recognized by observing the change in θ_d3. The θ_d3 increases quickly from θ₀₃ in the first stage, and keeps fluctuating around a quasi-steady level thereafter. The quasi-steady dynamics of fluids is a natural result of the quasi-constant viscous resistance, since the sum of FA and FB lengths in capillary remain unchanged and the inertial resistance is negligible after the very beginning stage. In fact, although the total fluid length in capillary is a constant, the difference in slip lengths between FA and FB should cause an increase in viscous resistance and a reduction in motion speed when FB is displaced by FA in Stage 2. Such a tendency turns out, however, not to have an obvious effect on the displacing curve in present simulation setup.

Fig. 3a demonstrates a typical result of spontaneous capillary displacement simulations for FA of different wetting properties (fixing A_B = −25, A_AB = −20, r = 10, while varying the A_A among −27.5, −30, −32.5, −35, −37.5, −40). Obviously, with the decrease of A_A (FA–capillary interaction enhances), the speed of the spontaneous capillary displacement process increases. As shown in Fig. 3b, the θ_d3 increases from the θ₀₃ and reaches a relatively high value, then fluctuates around a quasi-steady level. In case of starting by a large θ₀₃, the θ_d3 will remain at the quasi-steady level from the beginning. Generally speaking, it is found that the better the wetting property of FA, the smaller value the θ_d3 will finally reach, which corresponds to a larger driving force for the spontaneous capillary displacement process.

Fig. 3 Displacing curves (a) and the change in three-phase dynamic contact angle (θ_d3) (b) of the spontaneous capillary displacement process. Interaction parameter A_A varied among −40, −37.5, −35, −32.5, −30, and −27.5, while A_B = −25, A_AB = −20, r = 10.

Fig. 4a displays the simulations for FB of different wetting properties (fixing A_A = −40, A_AB = −20, r = 10, while varying A_B among −25, −27.5, −30, −32.5, −35, −37.5), and the variation in θ_d3 is shown in Fig. 4b. It is evident that, upon enhancing the wetting property of FB, the speed of the spontaneous capillary displacement process decreases and the quasi-steady level of θ_d3 increases. Such a tendency is as understandable as that of Fig. 3.

Fig. 4 Displacing curves (a) and the θ_d3 change (b) of the spontaneous capillary displacement process. Interaction parameter A_B varied among −37.5, −35, −32.5, −30, −27.5, and −25, which A_A = −40, A_AB = −20, r = 10.

Fig. 5a is the result for different FA–FB interaction (fixing A_A = −40, A_B = −30, r = 10, while varying A_AB among −15, −17.5, −20, −22.5, −25, −27.5), and the corresponding variation in θ_d3 is depicted in Fig. 5b. It is interesting to see that, upon increasing the miscibility between FA and FB in a certain broad range, the speed of the spontaneous capillary displacement process does not change significantly and the quasi-steady level of θ_d3 takes on a slow degradation.

Fig. 5 Displacing curves (a) and the θ_d3 change (b) of the spontaneous capillary displacement process. Interaction parameter A_AB varied among −27.5, −25, −22.5, −20, −17.5, and −15 while A_A = −35, A_B = −25, r = 10.

Fig. 6a shows the spontaneous capillary displacement results for different capillary radii (fixing A_A = −40, A_B = −30, A_AB = −20, while varying r among 5, 10, 15, 20, 25). It can be seen that, with the increase in capillary radius, the speed of the spontaneous capillary displacement process increases, so does the quasi-steady level of θ_d3 (Fig. 6b). Besides, the duration of initial stage of spontaneous capillary displacement seems to be related to the capillary radius; the greater the capillary radius, the longer the initial stage for the θ_d3 to reach its quasi-steady state.

Fig. 6 Displacing curves (a) and the θ_d3 change (b) of the spontaneous capillary displacement process. Capillary radius r varied among 5, 10, 15, 20, and 25, while A_A = −40, A_O = −30, A_AB = −20.

Consistency between spontaneous capillary displacement simulation and theoretical description

In the following, we verify the consistency between the above MDPD simulations and the theoretical description using eqn (9) for the spontaneous capillary displacement process.

We find that the simulation results can all match well the numerical calculations, provided appropriate parameters are employed in both simulation and calculation, among which the initial three-phase contact angle θ_i and the fluid front point are crucial.

The choice of the initial three-phase contact angle θ_i is important to the dynamics of the spontaneous capillary displacement process. An inappropriate beginning state will cause extra tensions on the meniscus and disturb the subsequent spontaneous capillary displacement process, resulting in an abnormal trace. To show the effect of θ_i, we have chosen four different values for comparison, i.e. 0°, 90°, 180° and θ₀₃. The choice of θ_i = 0° and 90° are popular among published studies, while θ_i = 180° is an imagined limiting case, where the interface is bended in an opposite direction. As displayed in Fig. 7a, only when θ_i is set to θ₀₃ can the simulation and the numerical calculation give consistent results; setting θ_i to 0°, 90°, or 180° will give deviated results. It seems the effect of θ_i on numerical calculation is greater than on MDPD simulation. Whereas numerical calculations give wide results when using inappropriate θ_i, the simulations produce parallel traces of fluid motion after an adjusting initial stage. This is understandable that, in MDPD simulations, if the fluid particles are not placed in equilibrium states, they are able to relax at the beginning of the spontaneous capillary displacement process.

Fig. 7 The consistency between MDPD simulation and numerical calculation using eqn (9) for the dynamic process of spontaneous capillary displacement. (a) The dependence of the three-phase initial contact angle θ_i. (b) The influence of the choice of fluid front.

One may be curious about the non-monotonic motion of the fluids in the case that θ_i equals 180°, which predicts the interface front first going in one direction and then the other. Note that, the arbitrarily formed initial meniscus renders a stress (could be calculated with Young–Laplace equation) that tend to propel the fluids to move in a reverse direction. Meanwhile, the difference in wetting property between FA and FB determines the driving force for the fluids to go in the forward direction. In some cases, this opposite stress is so great that, at the very beginning, the fluids will flow in a reverse direction, and after a relaxation for the bended meniscus, the flow would recover from this abnormal motion.

The choice of the FA–FB interface front is also a key parameter for the MDPD simulation to give results consistent with the numerical calculation. As illustrated in Fig. 7b, whereas taking the outmost or the central point, or somewhere in between, of the meniscus as the fluid front does not, in essence, affect the numerical calculation of the dynamics of spontaneous capillary displacement (the choice of fluid front is actually not a parameter in eqn (9)), it is important to the MDPD simulation. Only when taking the volume front as the front point of the meniscus, as we proposed in our previous article,³⁶ can the simulation results match completely the theoretical description.

Present numerical simulations have clearly show that multi-component fluid motions, such as fluid–fluid spontaneous capillary displacement within the fine capillary, can be well described by formula derived based on single-component fluid motion theory (i.e., Lucas–Washburn equation), which is consist with our recent experiment study¹² on both 2-phase (water, oil) and 3-phase (water, oil, gas) displacement. However it is worthy of noticing that, small but non-trivial modifications of the original Lucas–Washburn equation, as well as the way to pinning the interface front, are indispensable for precise tracking the fluid motions, although it remains challenging for experiments to capture the fluid movement at the very beginning.

Conclusion

In this article, we have conducted MDPD simulations on the spontaneous capillary displacement process. Various parameters, such as A_A, A_B, A_AB, and r, are adjusted separately to investigate their influence on the spontaneous capillary displacement. It is found that the displacing speed increases upon enhancing the FA wetting property, weakening the FB wetting property, or increasing the capillary radius; while the fluid–fluid miscibility seems not to be a sensitive factor in a certain range. We also find that the determination of the three-phase contact angle and the selection of the fluid front are crucial for the simulation, such that the simulated dynamic process of spontaneous capillary displacement can match, to a precise degree, the theoretical description based on numerical expression for fluid–fluid spontaneous capillary displacement (eqn (9)). We expect these findings would help to lay a basis for further simulation works on more complex systems, and bring us closer to practical experiments, or even to the development of rational strategies in industrial applications.

Acknowledgements

This work was financially supported by the National Basic Research Program (2012CB215503, 2012CB932800), the National Science Foundation of China (20933004, 21125312, 21303123), the National Hi-Tech R&D Program (2011AA050705), the Doctoral Fund of Ministry of Education of China (20110141130002), the Program for Changjiang Scholars and Innovative Research Team in University (IRT1030), and the PetroChina Changqing Oilfield Co.

References

1. J. F. Pankow and J. A. Cherry, Dense Chlorinated Solvents and Other DNAPLs in Groundwater: History, Behavior, and Remediation, Waterloo Press, Portland, Oregon, 1996.
2. J. W. Mercer and R. M. Cohen, A Review of Immiscible Fluids in the Subsurface: Properties, Models, Characterization and Remediation, J. Contam. Hydrol., 1990, 6, 107–163.
3. L. W. Lake, Enhanced Oil Recovery, Prentice Hall, Englewood Cliffs, New Jersey, 1989.
4. M. G. Gerritsen and L. J. Durlofsky, Modeling Fluid Flow in Oil Reservoirs, Annu. Rev. Fluid Mech., 2005, 37, 211–238.
5. H. Dehghanpour, B. Aminzadeh, M. Mirzaei and D. A. DiCarlo, Flow Coupling during Three-Phase Gravity Drainage, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 83, 065302(R).
6. M. I. J. Van Dijke and K. S. Sorbie, Pore-Scale Network Model for Three-Phase Flow in Mixed-Wet Porous Media, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2002, 66, 046302.
7. D. H. Fenwick and M. J. Blunt, Three-Dimensional Modeling of Three Phase Imbibition and Drainage, Adv. Water Resour., 1998, 21, 121–143.
8. T. Firincioglu, M. J. Blunt and D. Zhou, Three-Phase Flow and Wettability Effects in Triangular Capillaries, Colloids Surf., A, 1999, 155, 259–276.
9. H. Dehghanpour, B. Aminzadeh and D. A. DiCarlo, Hydraulic Conductance and Viscous Coupling of Three-Phase Layers in Angular Capillaries, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 83, 066320.
10. M. I. J. Van Dijke and K. S. Sorbie, The Relation between Interfacial Tensions and Wettability in Three-Phase Systems: Consequences for Pore Occupancy and Relative Permeability, J. Pet. Sci. Eng., 2002, 33, 39–48.
11. X. Xie, N. R. Morrow and J. S. Buckley, Contact Angle Hysteresis and the Stability of Wetting Changes Induced by Adsorption from Crude Oil, J. Pet. Sci. Eng., 2002, 33, 147–159.
12. W. Zhou, Y. Lu, C. Gao, W. Li, Y. Zhang, X. Li, C. Chen and J. Dong, Effects of Flow Pattern and Pore Size on Immiscible Continuous Three-Phase Displacement, Energy Fuels, 2013, 27, 717–724.
13. M. Sahimi, Flow Phenomena in Rocks: From Continuum Models to Fractals, Percolation, Cellular Automata, and Simulated Annealing, Rev. Mod. Phys., 1993, 65, 1393–1534.
14. K. S. Lee, N. Ivanova, V. M. Starov, N. Hilal and V. Dutschk, Kinetics of Wetting and Spreading by Aqueous Surfactant Solutions, Adv. Colloid Interface Sci., 2008, 144, 54–65.
15. N. V. Churaev, A. P. Ershov and Z. M. Zorin, Effect of Surfactants on the Kinetics of an Immiscible Displacement in Very Thin Capillaries, J. Colloid Interface Sci., 1996, 177, 589–601.
16. N. V. Churaev, A. P. Ershov, N. E. Esipova, R. M. Hill, V. D. Sobolev and Z. M. Zorin, Application of a Trisiloxane Surfactant for Removal of Oils from Hydrophobic Surfaces, Langmuir, 2001, 17, 1349–1356.
17. F. Tiberg, B. Zhmud, K. Hallstensson and M. Von Bahr, Capillary Rise of Surfactant Solutions, Phys. Chem. Chem. Phys., 2000, 2, 5189–5196.
18. C. D. Bain, Penetration of Surfactant Solutions into Hydrophobic Capillaries, Phys. Chem. Chem. Phys., 2005, 7, 3048–3051.
19. D. M. Colegate and C. D. Bain, Adsorption Kinetics in Micellar Solutions of Nonionic Surfactants, Phys. Rev. Lett., 2005, 95, 198302.
20. S. Girardo, R. Cingolani, S. Chibbaro, F. Diotallevi, S. Succi and D. Pisignano, Corner Liquid Imbibition during Capillary Penetration in Lithographically Made Microchannels, Appl. Phys. Lett., 2009, 94, 171901.
21. S. Girardo, S. Palpacelli, A. De Maio, R. Cingolani, S. Succi and D. Pisignano, Interplay between Shape and Roughness in Early-Stage Microcapillary Imbibition, Langmuir, 2012, 28, 2596–2603.
22. B. V. Zhmud, F. Tiberg. and J. Kizling, Dynamic Surface Tension in Concentrated Solutions of C_nE_m Surfactants: A Comparison between the Theory and Experiment, Langmuir, 2000, 16, 2557–2565.
23. P. S. Hammond and E. Unsal, Spontaneous and Forced Imbibition of Aqueous Wettability Altering Surfactant Solution into an Initially Oil-Wet Capillary, Langmuir, 2009, 25, 12591–12603.
24. P. S. Hammond and E. Unsal, Forced and Spontaneous Imbibition of Surfactant Solution into an Oil-Wet Capillary: The Effects of Surfactant Diffusion Ahead of the Advancing Meniscus, Langmuir, 2010, 26, 6206–6221.
25. P. S. Hammond and E. Unsal, Spontaneous Imbibition of Surfactant Solution into an Oil-Wet Capillary: Wettability Restoration by Surfactant-Contaminant Complexation, Langmuir, 2011, 27, 4412–4429.
26. A. Ghoufi, J. Emile and P. Malfreyt, Recent Advances in Many Body Dissipative Particles Dynamics Simulations of Liquid–Vapor Interfaces, Eur. Phys. J. E: Soft Matter Biol. Phys., 2013, 36, 10.
27. C. Chen, L. Zhuang, X. Li, J. Dong and J. Lu, A Many-Body Dissipative Particle Dynamics Study of Forced Water–Oil Displacement in Capillary, Langmuir, 2012, 28, 1330–1336.
28. W. L. Olbricht, Pore-Scale Prototypes of Multiphase Flow in Porous Media, Annu. Rev. Fluid Mech., 1996, 27, 187–213.
29. C. Cupelli, B. Henrich, T. Glatzel, R. Zengerle, M. Moseler and M. Santer, Dynamic Capillary Wetting Studied with Dissipative Particle Dynamics, New J. Phys., 2008, 10, 043009.
30. C. Chen, C. Gao, L. Zhuang, X. Li, P. Wu, J. Dong and J. Lu, A Many-Body Dissipative Particle Dynamics Study of Spontaneous Capillary Imbibition and Drainage, Langmuir, 2010, 26, 9533–9538.
31. R. Lucas, Ueber Das Zeitgesetz Des Kapillaren Aufstiegs Von Flussigkeiten, Kolloid-Z., 1918, 23, 15–22.
32. E. W. Washburn, The Dynamics of Capillary Flow, Phys. Rev., 1921, 17, 273–283.
33. P. Atkins and J. De Paula, Atkins' Physical Chemistry, Oxford University Press, Oxford, 8th edn, 2006.
34. P. De Gennes and G. Wetting, Statics and Dynamics, Rev. Mod. Phys., 1985, 57, 827–863.
35. S. P. Sutera and R. Skalak, The History of Poiseuille's Law, Annu. Rev. Fluid Mech., 1993, 25, 1–20.
36. C. Chen, K. Lu, L. Zhuang, X. Li, J. Dong and J. Lu, Effective Fluid Front of the Moving Meniscus in Capillary, Langmuir, 2013, 29, 3269–3273.
37. G. Martic, F. Gentner, D. Seveno, D. Coulon, J. De Coninck and T. D. Blake, A Molecular Dynamics Simulation of Capillary Imbibition, Langmuir, 2002, 18, 7971–7976.
38. M. R. Stukan, P. Ligneul, J. P. Crawshaw and E. S. Boek, Spontaneous Imbibition in Nanopores of Different Roughness and Wettability, Langmuir, 2010, 26, 13342–13352.
39. S. Gruener, T. Hofmann, D. Wallacher, A. V. Kityk and P. Huber, Capillary Rise of Water in Hydrophilic Nanopores, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 79, 067301.
40. W. Stroberg, S. Keten and W. K. Liu, Hydrodynamics of Capillary Imbibition under Nanoconfinement, Langmuir, 2012, 28, 14488–14495.
41. D. I. Dimitrov, A. Milchev and K. Binder, Capillary Rise in Nanopores: Molecular Dynamics Evidence for the Lucas–Washburn Equation, Phys. Rev. Lett., 2007, 99, 054501.
42. S. Supple and N. Quirke, Rapid Imbibition of Fluids in Carbon Nanotubes, Phys. Rev. Lett., 2003, 90, 214501.
43. S. Supple and N. Quirke, Molecular Dynamics of Transient Oil Flows in Nanopores I: Imbibition Speeds for Single Wall Carbon Nanotubes, J. Chem. Phys., 2004, 121, 8571–8579.
44. P. J. Hoogerbrugge and J. M. V. A. Koelman, Simulating Microscopic Hydrodynamic Phenomena with Dissipative Particle Dynamics, Europhys. Lett., 1992, 19, 155–160.
45. R. D. Groot and P. B. Warren, Dissipative Particle Dynamics: Bridging the Gap between Atomistic and Mesoscopic Simulation, J. Chem. Phys., 1997, 107, 4423–4435.
46. P. Español and P. Warren, Statistical Mechanics of Dissipative Particle Dynamics, Europhys. Lett., 1995, 30, 191–196.
47. A. A. Louis, P. G. Bolhuis and J. P. Hansen, Mean-Field Fluid Behavior of the Gaussian Core Model, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 2000, 62, 7961–7972.
48. P. B. Warren, Vapor–Liquid Coexistence in Many-Body Dissipative Particle Dynamics, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2003, 68, 066702.
49. J. L. Jones, M. Lal, J. N. Ruddock and N. A. Spenley, Dynamics of a Drop at a Liquid/Solid Interface in Simple Shear Fields: A Mesoscopic Simulation Study, Faraday Discuss., 1999, 112, 129–142.
50. J. A. Backer, C. P. Lowe, H. C. J. Hoefsloot and P. D. Iedema, Poiseuille Flow to Measure the Viscosity of Particle Model Fluids, J. Chem. Phys., 2005, 122, 154503.
51. S. Plimpton, Fast Parallel Algorithms for Short-Range Molecular Dynamics, J. Comput. Phys., 1995, 117, 1–19.
52. W. Humphrey, A. Dalke and K. Schulten, VMD – Visual Molecular Dynamics, J. Mol. Graphics, 1996, 14, 33–38.

---

Footnote

† Electronic supplementary information (ESI) available: Computational details and Tables S1 & S2. See DOI: 10.1039/c3ra47275b

---