A Many-Body Dissipative Particle Dynamics Study of Spontaneous

Español and Warren(21) showed that the dissipative and random forces should be related according ...... Prinsen , P., Warren , P. B., and Michels , M...
0 downloads 0 Views 2MB Size
pubs.acs.org/Langmuir © 2010 American Chemical Society

A Many-Body Dissipative Particle Dynamics Study of Spontaneous Capillary Imbibition and Drainage Chen Chen,† Chunning Gao,‡ Lin Zhuang,*,† Xuefeng Li,† Pingcang Wu,‡ Jinfeng Dong,† and Juntao Lu† † College of Chemistry and Molecular Sciences, Wuhan University, Wuhan 430072, China, and ‡Research Institute of Exploration and Development, PetroChina Changqing Oilfield Company, Xi’an 710018, China

Received January 8, 2010. Revised Manuscript Received February 16, 2010 The spontaneous capillary imbibition and drainage processes are studied using many-body dissipative particle dynamics (MDPD) simulations. By adjusting the solid-liquid interaction parameter, different wetting behavior between the fluid and the capillary wall, corresponding to the static contact angle ranging from 0° to 180°, can be controllably simulated. For wetting fluids, the spontaneous capillary imbibition (SCI) is evident in MDPD simulations. It is found that, whereas the corrected Lucas-Washburn equation (taking into account the dynamic contact angle and the fluid inertia) can well describe the SCI simulation result for the completely wetting fluid, it deviates, to a notable degree, from the results of partly wetting fluids. In particular, this corrected equation cannot be used to describe the spontaneous capillary drainage (SCD) processes. To solve this problem, we derive an improved form of the Lucas-Washburn equation, in which the slip effects of fluid particles at the capillary wall are treated. Such an improved equation turns out to be capable of describing all the simulation results of both the SCI and the SCD. These findings provide new insights into the SCI and SCD processes and improve the mathematical base.

Introduction The spontaneous capillary imbibition (SCI) is a ubiquitous phenomenon in nature.1,2 Not only can SCI be seen in many biological processes, such as the sap driving in plants3 and transportations across biomembranes,4 but it has also been widely applied in modern technologies, such as oil recovery,5 ink printing,6 dyeing of textiles,7 powder dissolution,8 and nanofluidic devices.9 Over almost a century, the SCI phenomena have been studied experimentally and computationally. According to experimental measurements, the SCI process was found to be well describable by the Lucas-Washburn equation:10,11 sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi γrt cos θ ð1Þ l ¼ 2η where l is the penetration length of the fluid, γ its vapor-liquid surface tension, and η its shear viscosity; r is the pore radius of the capillary; θ is the equilibrium contact angle between the fluid and the capillary wall; and t is the time elapsed. Yet most of the SCI experiments were conducted at macroscopic scales;12,13 for systems of much smaller scales, where *Corresponding author. E-mail: [email protected].

(1) Alava, M.; Dube, M.; Rost, M. Adv. Phys. 2004, 53, 83. (2) De Gennes, P. G. Rev. Mod. Phys. 1985, 57, 827. (3) Pickard, W. F. Prog. Biophys. Mol. Biol. 1982, 37, 181. (4) Bramlage, W. J.; Leopold, A. C.; Parrish, D. J. Plant Physiol. 1978, 61, 525. (5) Morrow, N. R.; Mason, G. Curr. Opin. Colloid Interface Sci. 2001, 6, 321. (6) Clarke, A.; Blake, T. D.; Carruthers, K.; Woodward, A. Langmuir 2002, 18, 2980. (7) Pezron, I.; Bourgain, G.; Quere, D. J. Colloid Interface Sci. 1995, 173, 319. (8) Brielles, N.; Chantraine, F.; Viana, M.; Chulia, D.; Branlard, P.; Rubinstenn, G.; Lequeux, F.; Lasseux, D.; Birot, M.; Roux, D.; Mondain-Monval, O. Ind. Eng. Chem. Res. 2007, 46, 5785. (9) Supple, S.; Quirke, N. Phys. Rev. Lett. 2003, 90, 214501. (10) Lucas, R. Kolloid Z. 1918, 23, 15. (11) Washburn, E. W. Phys. Rev. 1921, 17, 273. (12) Zhmud, B. V.; Tiberg, F.; Hallstensson, K. J. Colloid Interface Sci. 2000, 228, 263. (13) Kornev, K. G.; Neimark, A. V. J. Colloid Interface Sci. 2003, 262, 253.

Langmuir 2010, 26(12), 9533–9538

experimental observations become difficult and computational simulations are playing a major role, the validity of the LucasWashburn equation has been challenged.9,14-16 Martic et al.14 had pointed out that, for microscopic systems, the contact angle θ in the equation should be taken as the dynamic contact angle θd rather than the static contact angle θ0, and the inertia contribution should also be taken into account. Although such a correction in the Lucas-Washburn equation turns out to be successful for completely wetting fluids, which is assumed to follow Poiseuille’s law17 and results in zero fluid speed at the capillary wall (the socalled “no-slip” boundary conditions), it would deviate from the situation of partially wetting or nonwetting fluids because in these cases the solid-liquid interactions decrease to a significant degree and the fluid velocity at the interface cannot be neglected. In the present work, we studied the SCI and the spontaneous capillary drainage (SCD) using many-body dissipative particle dynamic (MDPD) simulations. We found that the LucasWashburn equation clearly deviate from the MDPD simulation results of either SCI or SCD processes, but by introducing a new parameter into the equation to correct the fluid velocity at the capillary wall, the improved Lucas-Washburn equation can describe very well both the SCI and SCD behavior of fluids with distinct wetting properties.

Computational Method There has been a robust growth in the development of simulation methods for fluids, among which is the dissipative particle dynamics (DPD) method, specially designed for fluid problems at micro- and mesoscales. The DPD was originally proposed by Hoogerbrugge and Koelman,18 with the purpose of simulating (14) Martic, G.; Gentner, F.; Seveno, D.; Coulon, D.; De Coninck, J. Langmuir 2002, 18, 7971. (15) Supple, S.; Quirke, N. J. Chem. Phys. 2004, 121, 8571. (16) Dimitrov, D. I.; Milchev, A.; Binder, K. Phys. Rev. Lett. 2007, 99, 054501. (17) Sutera, S. P.; Skalak, R. Annu. Rev. Fluid Mech. 1993, 25, 1. (18) Hoogerbrugge, P. J.; Koelman, J. M. V. A. Europhys. Lett. 1992, 19, 155.

Published on Web 03/12/2010

DOI: 10.1021/la100105f

9533

Article

Chen et al.

complex fluids at a coarse-grained level, and then further developed by several investigators in its physical essence.19 The DPD method can be considered as a molecular dynamics (MD) method with two important innovations:20 using relatively soft interaction potentials and employing a locally momentum conserved thermostat. Particles in DPD are placed in a continuous space and characterized by its position ri, velocity vi, and mass mi, and the movements are governed by Newton’s second law. Forces acting on each particle are separated into the conservative, dissipative, and random parts: R Fij ¼ FCij þ FD ij þ Fij

ð2Þ

and are designed to vanish beyond a certain cutoff distance rc (in the reduced unit of length). The conservative part is soft repulsive force acting along the center of particles: FCij ¼ AwC ðrij Þeij

ð3Þ

where A is the maximum repulsion between particle i and j, wC(rij) is a weight function vanishing for r > rc = 1, and rij = rj - ri, rij = |rij|, and eij = rij/rij. The dissipative part and random part are given by FD ij ¼ -γwD ðrij Þðvij 3 eij Þeij FR ij

¼ σwR ðrij Þξij eij

FCij ¼ Aij wC ðrij Þeij þ BðF i þ F j Þwd ðrij Þeij

ð6Þ

in which the originally repulsive soft potential is tuned to be attractive by setting A < 0 and a local density dependent manybody repulsive force is added by setting B > 0, while F is the partial amplitude depending on a weighted local density which is defined for the ith particle to be Fi ¼

X

wF ðrij Þ

ð7Þ

j6¼i

wF ðrÞ ¼ 15ð1 -r=rd Þ2 =ð2πrd 3 Þ

ð8Þ

The two weight functions wC and wd in eq 6 are in similar definitions: wC ðrÞ ¼ 1 -r=rc

ð9Þ

wd ðrÞ ¼ 1 -r=rd

ð10Þ

ð4Þ ð5Þ

where γ is a coefficient controlling the extent of dissipation, ξij is a random number drawn from a Guassian distribution of zero mean and variance equal to 1/δt, and wD(rij) and wR(rij) are weight functions vanishing for r > rc = 1. Espa~nol and Warren21 showed that the dissipative and random forces should be related according to wD(rij) = [wC(rij)]2 and σ2 = 2γkBT, such that the dissipative and random forces behave as a thermostat without influencing the equilibrium thermodynamics of the system. Although DPD performed excellent in many fluid-related problems, such as copolymer melts,22,23 surfactant solutions,24,25 phase separation kinetics,26,27 and the fusion of biomembranes,28 it is well recognized that the quadratic equation of state (EOS) produced by the conservative force is quite different from the EOS of real fluids which contains a van der Waals loop.29 Simply speaking, the solely repulsive interaction between particles drives the DPD beads to fill the whole simulation box, making it impossible to study systems containing free surfaces, such as vapor-liquid interfaces and fluid flows in capillaries. Until recently, several groups20,30,31 have made efforts in modifying the original DPD method to include the van der Waals loop in the EOS so as to make the DPD suitable for simulations of fluid systems with free surfaces. These augmented DPD methods are (19) Groot, R.; Warren, P. B. J. Chem. Phys. 1997, 107, 4423. (20) Warren, P. B. Phys. Rev. E 2003, 68, 066702. (21) Espa~nol, P.; Warren, P. Europhys. Lett. 1995, 30, 191. (22) Groot, R. D.; Madden, T. J.; Tildesley, D. J. J. Chem. Phys. 1999, 110, 9739. (23) Groot, R. D. Langmuir 2000, 16, 7493. (24) Jury, S.; Bladon, P.; Cates, M.; Krishna, S.; Hagen, M.; Ruddock, N.; Warren, P. Phys. Chem. Chem. Phys. 1999, 1, 2051. (25) Prinsen, P.; Warren, P. B.; Michels, M. A. J. Phys. Rev. Lett. 2002, 89, 148302. (26) Groot, R. D.; Madden, T. J. J. Chem. Phys. 1998, 108, 8713. (27) Warren, P. B. Phys. Rev. Lett. 2001, 26, 225702. (28) Shillcock, J. C.; Lipowsky, R. Nat. Mater. 2005, 4, 225. (29) Louis, A. A.; Bolhuis, P. G.; Hansen, J. P. Phys. Rev. E 2000, 62, 7961. (30) Pagonabarraga, I.; Frenkel, D. J. Chem. Phys. 2001, 115, 5015. (31) Trofimov, S. Y.; Nies, E. L. F.; Michels, M. A. J. J. Chem. Phys. 2002, 117, 9383.

9534 DOI: 10.1021/la100105f

called MDPD, with M denoting many-body or multibody. In the present work, we employ the MDPD method reported by Warren,20 in which both the random and dissipative parts of the force are kept unchanged while conservative part is represented as

but with different cutoff distances: rc = 1 and rd = 0.75. With these modifications, both planar and curved vapor-liquid interfaces can be achieved and in excellent agreement with the theory.20 In the present work, the attraction parameters for solid-solid and liquid-liquid interactions, Ass and All, are fixed as -40.0, while leaving the solid-liquid interaction Asl tunable. All the repulsion parameters B are kept as 25, the amplitude of noise σ is set to be 6.0, and the time step in all simulations was set to be 0.01 MDPD time units. Under such conditions, the equilibrium number density F for fluid is found to be about 6.0. We also calculated the surface tension and viscosity following the methods described in detail elsewhere,32,33 which give the value γ = 7.62 ( 0.05 for a 10  10  10 slab containing 6083 particles in a 20  10  10 box and η = 7.39 ( 0.06 for a 12  8  8 box filled with 4604 particles. Both simulation data were obtained with first 1000 MDPD time units for equilibrium and then 4000 MDPD time units for data production. In summary, the wetting behavior of the fluid in the capillary can be tuned by only adjusting one variable, the solid-liquid interaction parameter Asl. In the present work, the MDPD algorithm was implemented by modifying the DPD code in the LAMMPS package.34

Results and Discussion Static State of Fluids in Capillary. Before studying the dynamics of fluids in the capillary, the static state should be first simulated so as to obtain the static contact angle, which defines the wetting property of the fluid on the capillary wall. The system we use for the calculation of static contact angles is similar to that reported by Cupelli et al.:33 a half-sealing capillary partially filled with a fluid (Figure 1). The system was built in two steps: First, in (32) Jones, J. L.; Lal, M.; Ruddock, J. N.; Spenley, N. A. Faraday Discuss. 1999, 112, 129. (33) Backer, J. A.; Lowe, C. P.; Hoefsloot, H. C. J.; Iedema, P. D. J. Chem. Phys. 2005, 122, 154503. (34) Plimpton, S. J. J. Comput. Phys. 1995, 117, 1.

Langmuir 2010, 26(12), 9533–9538

Chen et al.

Article

Figure 2. Relationship between the solid-liquid interaction para-

Figure 1. Simulation system for static contact angle calculations. (a) Schematic illustration of dividing the fluid into concentric cylinders with certain thickness. (b) and (c) are the slide view of the capillary at the initial and final state of simulation, respectively. Brown balls are the slab particles, black balls the capillary particles, and gray balls the fluid particles.

meter Asl and the equilibrium static contact angle θ0 obtained by MDPD simulations.

a 40  12  12 simulation box a group of 34 560 homogeneous particles were arranged in a capillary shape and relaxed for 100 time units to equilibrate the particle positions, then the particles at the left-hand end of the capillary were defined as the sealing slab, the particles at the outer surface of the capillary were defined as the capillary wall, and the rest particles (in the capillary) are the fluid. Second, a half of the fluid particles at the right-hand side were deleted so as to leave space for the development of the meniscus fluid surface (Figure 1b). We use the model described by Henrich et al.35 to treat the particles of the capillary wall and the sealing slab, i.e., pinning them at their equilibrated positions by a spring force: Fs ¼ ks ðri ðtÞ -ri ð0ÞÞ

ð11Þ

with the elastic coefficient ks = 12.5. This model was proved to be able to produce a solid-liquid interface without considerable oscillations in the density and temperature, in comparison to those rigid wall models.35 Simulations were then further taken in a NVE ensemble, and the static contact angle was calculated after equilibrating the thusprepared system for 1000 time units (Figure 1c). To extract the contact angle, we divided the fluids into concentric cylindrical shells with certain thickness (Figure 1a, Δr = 1rc in this work) and then exported the x position of the front for each shell; by fitting all the x position to a cyclic function of z, the contact angle resulted. Considering its special role (vide infra), the outmost shell in contact with the capillary wall was omitted in data fitting.14,16 As shown in Figure 2, by tuning the solid-liquid interaction parameter Asl from -40 to -2.5, different static contact angles θ0 were obtained, ranging from 0° to 180°. In other words, it is controllable to simulate the wetting behavior of a fluid in a capillary, from completely solvophilic to completely solvophobic. In the following simulations of SCI and SCD, such a θ0-Asl relationship was used to define the wetting property of a fluid in the capillary. SCI and SCD of Completely Wetting Fluid. The SCI and SCD processes were simulated by putting a capillary array (with each in radius of 10 and consisting of 39 438 wall particles) onto a (35) Henrich, B.; Cupelli, C.; Moseler, M.; Santer, M. Europhys. Lett. 2007, 80, 60004.

Langmuir 2010, 26(12), 9533–9538

Figure 3. Simulation system for the imbibition and drainage processes. (a) and (b) are the top view and the side view, respectively, of the capillary array (periodic in y and z directions) putting onto a fluid reservoir (in gray), lr and lc are fluid lengths in the reservoir and the capillary, respectively. (c) and (d) are the section view of the capillary at the initial and final state of SCI simulations, respectively. (e) and (f) are the section view of the capillary at the initial and final state of SCD simulations, respectively.

fluid reservoir (a thick liquid layer of 112 491 particles, 38  22  22 in x, y, and z dimensions for SCI, while 58 273 particles, 20  22  22 in x, y, and z dimensions for SCD). As shown in Figure 3a,b, the system was designed to be periodic in the y and z directions, and the imbibition/drainage process occurred in the capillaries at the x direction. In the simulation of SCI, the initial meniscus position was set at the entrance of the capillary with a static contact angle (Figure 3c), and the imbibition process was recorded afterward (Figure 3d), while the SCD process started from a fully filled state (Figure 3e) and was recorded during the subsequent receding process (Figure 3f). Upon tuning the interaction between the fluid and the capillary wall (by setting the Asl parameter), the imbibition/drainage process occurred at different rates, characterized by the shift in the meniscus after 400 time units of simulation. As shown in Figure 4, in the range of Asl = -27.5 to -40 (corresponding to θ0 changing from slightly below 90° to 0°), the SCI process can be observed; when Asl > -26.5 (θ0 > 90°), the SCD process occurs. The DOI: 10.1021/la100105f

9535

Article

Chen et al.

Figure 4. Relationship between the solid-liquid interaction parameter Asl and the shifting length of the fluid meniscus after 400 time units of MDPD simulation. The left part is the SCI result and right part the SCD result.

Figure 5. Development of the fluid meniscus, as a function of time, in the SCI simulation for a completely wetting fluid (Asl = -40). Theoretical predictions of the original and corrected LucasWashburn equations are also plotted for comparison.

shifting rate of the fluid meniscus increases, as expected, upon enhancing the solvophilicity or solvophobicity. The dynamics of SCI has been extensively studied by experiments and can be well described by the Washburn equation: P dl r2 P ¼ 8ηl dt

ð12Þ

where dl/dt is the shifting rate of the fluid in the capillary, l the shifting distance, t the observation time, P r the radius of the capillary, η the viscosity of fluid, and P the total effective pressure driving the fluid along the capillary. In the case of spontaneous imbibition in horizontal capillaries, this formula can be reduced to eq 1, the Lucas-Washburn equation. As mentioned in the Introduction, Cupelli et al.36 suggested to correct the Lucas-Washburn equation by taking into account the dynamic contact angle θd and the inertia, which leads to the following equation (assuming the capillary to be in a round shape): 1 dpðtÞ 2γlv 8η dl ¼ cosðθd Þ - 2 l πr2 dt r dt r

ð13Þ

where dp(t) is the momentum change of fluid during time interval dt, γlv the liquid-vapor surface tension of fluid, l the fluid length within capillary in the x direction. Further, the dp(t)/dt term was separated into two parts:   dpðtÞ d dlc dlr ¼ mc ðtÞ þ Rmr ðtÞ dt dt dt dt

ð14Þ

where m is the mass of fluid, with the c and r in the subscript denoting the capillary and the fluid reservoir, respectively. In addition to the momentum change of the fluid in the capillary, the momentum change of the fluid reservoir (the left part in Figure 3b) was also considered and weighted with a parameter R. While the physical meaning of R is not straightforward at the moment and remains a topic for further study, we find that incorporating R is necessary for acceptable data fitting. Figure 5 presents the MDPD simulation of the SCI process of a completely wetting fluid (Asl = -40), along with the predictions (36) Cupelli, C.; Henrich, B.; Glatzel, T.; Zengerle, R.; Moseler, M.; Santer, M. New J. Phys. 2008, 10, 043009.

9536 DOI: 10.1021/la100105f

Figure 6. MDPD simulation results of the imbibition process of a completely wetting fluid (Asl = -40) with different thickness l0 (in MDPD length unit) of the fluid reservoir. Inset: the initial stage of SCI processes.

of Lucas-Washburn equation and eq 13 with different R values. Clearly, the MDPD simulation result cannot be fitted by the original Lucas-Washburn equation, but by using dynamic contact angles instead of the static contact angle, the MDPD result can be described except for the initial stage of SCI process. Only when the momentum change in the fluid reservoir is also taken into account (R = 0.5 turns out to be the best in this work) can eq 13 closely match the full MDPD result. To further reveal the dependence of the initial stage of SCI process on the momentum of the fluid reservoir, we carried out the above MDPD simulation with different volumes of the fluid reservoir. As demonstrated in Figure 6, due to the inertia difference, the imbibition rates are not the same in the initial stage but converge gradually. Except for the difference in the initial rate, all of the SCI processes can be well fitted by eq 13 with the same R value, i.e., R = 0.5. It seems R is mostly determined by the area fraction of the capillary at the capillary/fluid reservoir interface. It was found, however, that neither the initial nor the final states of the SCD process were dependent on R. Slip Effects in SCI and SCD Processes. Although eq 13 turns out to be capable of describing the MDPD result of the completely wetting fluid, it deviates for the other cases where Asl is not -40, i.e., θ0 is not 0° (vide infra). Such an observation should be ascribed to the inadequate assumption of the “no-slip” Langmuir 2010, 26(12), 9533–9538

Chen et al.

Article

Figure 7. Schematic illustration of the definition of slip velocity and slip length. The fluid retains a nonzero slip velocity (vs in red line) at the solid-liquid boundary, which can also be measured in terms of slip length b, a virtual depth obtained by extrapolating the velocity profile into the solid to meet the zero-velocity boundary condition.

boundary condition (zero velocity of fluid particles at the capillary wall) when applying Poiseuille’s law to derive the LucasWashburn equation. This assumption should, in principle, hold for completely wetting fluids; i.e., the interaction between the fluid and the capillary wall is extremely strong, but it may not be the case for partially wetting or nonwetting fluids, where the slip effects at the fluid/capillary interface cannot be neglected. Dimitrov et al.16 have suggested replacing the capillary radius in the Lucas-Washburn equation with an effective radius so as to involve the slip effects. Here we try to derive the equation more strictly in mathematics. In Navier’s model,37 the slip effect is characterized by the slip length b or the slip velocity vs, the definitions of which are illustrated in Figure 7. The b and vs are interchangeable through  Dv vs ¼ -b  Dz

ð15Þ boundary

The slip length b can be used to describe the deviation of practical boundary conditions from the “no-slip” assumption. Taking the slip effect into consideration, eq 13 can be modified as 1 dpðtÞ 2γlv 8η dl ¼ cosðθd Þ - 2 l πr2 dt r þ 4br dt r

ð16Þ

The MDPD simulations for the SCI of partially wetting fluids and the SCD of nonwetting fluids are exemplified in Figure 8, along with the curve fitting using eqs 13 and 16. It is evident that, while the slip effect plays a subtle role in the SCI process, it becomes significant in the SCD. This observation is understandable because the interaction between the nonwetting fluid and the capillary wall is so small that the fluid particles at the capillary wall have to move at a velocity comparable to that of the bulk particles. It is also imaginable that the solid-liquid interaction parameter Asl and the slip length b should be related; the weaker the interaction, the greater the slip length. Such a deduction is verified by the b-Asl relationship shown in Figure 9, which was obtained by simulating the SCI and SCD processes of fluids with different Asl and then fitting the data to eq 16. We did not get reliable results for Asl ranging from -32.5 to -25 because of the fluctuation in the dynamic contact angle around 90°. Nevertheless, the b-Asl relationship seems to be describable by a quadratic curve (the red dashed line in Figure 9). The slip effect also has an interesting impact on the dynamic contact angle in SCI and SCD. As depicted in Figure 10, for (37) Navier, C. L. M. H. Mem. Acad. Sci. Inst. Fr. 1827, 6, 839.

Langmuir 2010, 26(12), 9533–9538

Figure 8. Fluid length in the capillary, changing with time, during the SCI simulation for a partly wetting fluid (Asl = -32.5) and the SCD simulation for a nonwetting fluid (Asl = -20). Theoretical predictions of the corrected Lucas-Washburn equation (eq 13) and our improved equation (eq 16) are also plotted for comparison.

Figure 9. Relationship between the solid-liquid interaction parameter Asl and the slip length b obtained by fitting the SCI and SCD simulation results to eq 16. This relationship seems to fit a quadratic curve (red dashed line).

Figure 10. Change in the dynamic contact angle θd during the SCI and SCD processes. Completely wetting (Asl = -40, red line), partly wetting (Asl = -35, black line), and nonwetting (Asl = -20, blue line) examples are presented.

completely and partially wetting fluids (e.g., Asl = -40, θ0 = 0° and Asl = -35, θ0 = 46.8°, respectively), the θd increases rapidly at the initial stage of SCI and then decreases slowly afterward, and the θd remains noticeably greater than the corresponding θ0 in the DOI: 10.1021/la100105f

9537

Article

period of MDPD simulations. This is because the fluid-capillary interaction causes a friction effect resisting the spreading of the fluid during the SCI process. In SCD, however, the fluid-capillary interaction is weak, and the θd was observed to fluctuate just around θ0 in the MDPD simulations.

Conclusion The spontaneous capillary imbibition (SCI) and spontaneous capillary drainage (SCD) were simulated in the present work using the many-body dissipative particle dynamic (MDPD) method. By adjusting the solid-liquid interaction parameter Asl, it is feasible to simulate the wetting behavior of a fluid in the capillary from completely wetting to completely nonwetting, corresponding to the change in the static contact angle θ0 ranging from 0° to 180°. On the basis of such a control in the wetting behavior, the SCI of wetting fluids and SCD of nonwetting fluids were studied. It was found that, whereas the MDPD simulation result of the SCI process of completely wetting fluid is in consistent with the prediction of the corrected Lucas-Washburn equation (eq 13),

9538 DOI: 10.1021/la100105f

Chen et al.

the MDPD results for partially wetting and nonwetting fluids cannot be well described by this known equation, and the slip effects of the fluid particles at the solid-liquid interface have to be taken into account. By incorporating the concept of slip length into the LucasWashburn equation to correct for the nonzero velocity of fluid particles at the solid-liquid interface, we have derived a generalized equation for fluids with any wetting property. Such an improved equation can describe, to an excellent degree, all the MDPD simulation results for both the SCI and the SCD processes, and new insights, such as the b-Asl relationship, have thus been gained. The particular capability of our equation for the study of SCD, a process as important as but much less investigated than the SCI, makes it superior to the Lucas-Washburn equation and would benefit studies in this field. Acknowledgment. This work was financially supported by the National Science Foundation of China (Grant 20933004) and the PetroChina Changqing Oilfield Co.

Langmuir 2010, 26(12), 9533–9538