Molecular dynamics of wetting layer formation and forced water invasion in angular nanopores with mixed wettability.

Molecular dynamics of wetting layer formation and forced water invasion in angular nanopores with mixed wettability Mohammad Sedghi, Mohammad Piri, and Lamia Goual Citation: The Journal of Chemical Physics 141, 194703 (2014); doi: 10.1063/1.4901752 View online: http://dx.doi.org/10.1063/1.4901752 View Table of Contents: http://scitation.aip.org/content/aip/journal/jcp/141/19?ver=pdfcov Published by the AIP Publishing Articles you may be interested in Molecular dynamics simulations of wetting behavior of water droplets on polytetrafluorethylene surfaces J. Chem. Phys. 140, 114704 (2014); 10.1063/1.4868641 Wetting kinetics of water nano-droplet containing non-surfactant nanoparticles: A molecular dynamics study Appl. Phys. Lett. 103, 253104 (2013); 10.1063/1.4837717 Molecular dynamics studies of water deposition on hematite surfaces AIP Conf. Proc. 1504, 780 (2012); 10.1063/1.4771809 Analysis on wetting and local dynamic properties of single water droplet on a polarized solid surface: A molecular dynamics study J. Chem. Phys. 135, 014703 (2011); 10.1063/1.3601055 Volumetric Fraction Dynamic Measurement in Oil‐Water‐Gas Multiphase Horizontal Pipe Flow with Dual Energy Gamma‐Ray AIP Conf. Proc. 914, 409 (2007); 10.1063/1.2747460

This article is copyrighted as indicated in the article. Reuse of AIP content is subject to the terms at: http://scitation.aip.org/termsconditions. Downloaded to IP: 130.179.16.201 On: Tue, 23 Jun 2015 06:58:27

THE JOURNAL OF CHEMICAL PHYSICS 141, 194703 (2014)

Molecular dynamics of wetting layer formation and forced water invasion in angular nanopores with mixed wettability Mohammad Sedghi,a) Mohammad Piri, and Lamia Goual Department of Chemical and Petroleum Engineering, University of Wyoming, 1000 East University Avenue, Laramie, Wyoming 82071, USA

(Received 11 August 2014; accepted 3 November 2014; published online 19 November 2014) The depletion of conventional hydrocarbon reservoirs has prompted the oil and gas industry to search for unconventional resources such as shale gas/oil reservoirs. In shale rocks, considerable amounts of hydrocarbon reside in nanoscale pore spaces. As a result, understanding the multiphase flow of wetting and non-wetting phases in nanopores is important to improve oil and gas recovery from these formations. This study was designed to investigate the threshold capillary pressure of oil and water displacements in a capillary dominated regime inside nanoscale pores using nonequilibrium molecular dynamics (NEMD) simulations. The pores have the same cross-sectional area and volume but different cross-sectional shapes. Oil and water particles were represented with a coarse grained model and the NEMD simulations were conducted by assigning external pressure on an impermeable piston. Threshold capillary pressures were determined for the drainage process (water replaced by oil) in different pores. The molecular dynamics results are in close agreements with calculations using the Mayer-Stowe-Princen (MS-P) method which has been developed on the premise of energy balance in thermodynamic equilibrium. After the drainage simulations, a change in wall particles’ wettability from water-wet to oil-wet was implemented based on the final configuration of oil and water inside the pore. Waterflooding simulations were then carried out at the threshold capillary pressure. The results show that the oil layer formed between water in the corner and in the center of the pore is not stable and collapses as the simulation continues. This is in line with the predictions from the MS-P method. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4901752] I. INTRODUCTION

Recent advancements in horizontal drilling and hydraulic fracturing technologies have made it possible to commercially produce oil and gas from shale and ultra-tight sandstone formations. The pore diameters in these rocks are in the nanometer range, which results in very low porosities and permeabilities. Considering nanoscale phenomena such as slip boundary conditions, non-continuum effects, higher ratio of wall surface area to pore volume (PV), etc., one would expect different fluid flow and phase behavior inside nanoscale confinements than in microscale and larger pores.1–3 Developing an improved understanding of the fundamentals of multiphase flow in nano-porous media is crucial to enhance oil recovery from shales; One of the key aspects of multiphase flow in nanopores that has not been investigated in the past is the effect of wall-fluids interactions on the threshold capillary pressure of fluid/fluid displacements in these systems. Most of the molecular dynamics (MD) studies on fluid flow in nanopores have considered nanotubes with circular cross sections.4–16 However, the majority of pores in porous media have sharp corners and could be best represented with angular cross sections such as triangles and squares.18, 19 The difference between polygonal and circular cross sections is that the wetting phase will retract to and remain in the corners as the non-wetting phase starts to occupy the center of the a) [email protected]

0021-9606/2014/141(19)/194703/12/$30.00

pore during a fluid-fluid displacement. For instance, in a displacement of water by oil in an initially water-wet pore during a drainage process, if the pore cross section is circular, then a thin film of water will remain on the pore wall. However, if the cross section is angular, the water will leave thick layers in the corners as well as a thin film on the mid-section of the pore wall similar to the case in circular pores. Under certain conditions, the water film may rupture leading to a direct contact between the oil and the pore walls. As a result, the wettability of portions of the wall might change from water-wet to weakly or strongly oil-wet after a sufficient amount of time. This is due to the interactions between heavier molecules of oil called resins and asphaltenes with polar components of the wall surface. This mixed wettability has been extensively characterized in the literature for macropores.17–20 After primary drainage, waterflooding can occur by displacing the oil, that occupies the center of the pore, by water. The configuration of water and oil at the end of the waterflooding process can vary depending on the wettability condition of the pore walls. If the walls remained uniformly waterwet, meaning that either the water film in the middle part of the pore did not rupture or it was impermeable to oil components, then water would occupy the entire pore space after waterflooding. On the other hand, if the mid-section of the pore becomes strongly oil-wet, then a thin layer of oil could form between the water in the corners and the water in the center of the pore after waterflooding. The formation of the wetting layer depends on the interfacial properties of the oil,

141, 194703-1

© 2014 AIP Publishing LLC


194703-2

Sedghi, Piri, and Goual

water, and solid surface as well as the initial saturation of the water and the pressure difference between the fluid phases. The dynamics of fluids flow inside nanopores starts to deviate from the continuum model equations when the pore dimension becomes comparable to the molecular size.3, 21 Hence, a practical tool to study flow and transport under these conditions is MD simulation. Over the last decade, the spontaneous imbibition of wetting phases into nanotubes has been extensively investigated through MD simulations. One of the early studies in this area was performed by Martic et al.6 who modified the Lucas-Washburn equation by replacing the equilibrium contact angle with the dynamic contact angle. In their equations, the velocity of the rising meniscus was derived from the molecular kinetic theory. They validated the modified Lucas-Washburn equation by using MD simulations of a Lennard-Jones (LJ) fluid spontaneously imbibing a nanotube. The position of the advancing meniscus and the dynamic contact angle were measured during the simulation and the results followed the trend predicted by the modified Lucas-Washburn equation. Their work was later followed by Stukan et al.4 who addressed the effects of roughness and wettability of the wall on the dynamics of spontaneous imbibition. In another study, Dimitrov et al.9 obtained the velocity of capillary rise under an external pressure in order to compute the permeability and the capillary pressure of the system. The external pressure was applied to an impermeable piston placed at the end of the simulation box. Contrary to the above-mentioned studies that confirmed the validity of Lucas-Washburn equation for spontaneous imbibition, MD studies performed by Supple and Quirke12 revealed that the velocity of decane meniscus rise in graphene nanotubes is proportional to the time instead of its square root. Considering their simulation times were less than 1 ns, it is possible that they had observed initial dynamics effects ignored in other studies. Perhaps, one of the most relevant studies on two-phase flow in nanopores is the one presented by Chen et al.15 where forced displacement of oil by water was investigated in a circular nanotube and the amount of remaining oil and the threshold external force (presented in reduced units) where determined for different pore wettabilities. The results indicated that the amount of residual oil and the threshold external force increase as the pore wall becomes more oil-wet. In this study, we focus on the displacements of oil and water inside nanopores with angular cross sections and with varying pore wall wettabilities. We investigate the threshold capillary pressures and pore fluid configurations associated with drainage and forced waterflooding displacement processes. Formation of wetting layers of water in the corners of water-wet pores, establishment of strongly oil-wet pore walls through wettability alteration process, and formation and collapse of wetting layers of oil in oil-wet pores are studied. To the best of our knowledge, this is the first time that molecular dynamics of wetting layer formation and forced water-oil displacements in angular nanopores with mixed wettabilities are examined. This paper is organized as follows: In Sec. II, we describe our systems of study and the MD parameters and techniques as well as the thermodynamic equations derived using the Mayer-Stowe-Princen (MS-P) method. In Sec. III, we compare the threshold capillary pressure obtained from

J. Chem. Phys. 141, 194703 (2014)

MD simulations with calculations from the MS-P method for both drainage and waterflooding processes. A set of final remarks are summarized in the Conclusions. II. METHODOLOGY A. Threshold capillary pressure

In this work, we used the definition of the capillary pressure as the pressure differential between adjacent oil and water phases inside capillary pores. With this definition, capillary pressure can be positive or negative depending on the pore wettability, although only the absolute values are reported in this paper. For a circular capillary tube with an inside radius of R, the threshold capillary pressure for the nonwetting phase j to invade the tube and displace the wetting phase i, is calculated from the Young-Laplace equation Pcthij = Pi − Pj =

2γij cos(θij ) R

,

(1)

where γ ij is the interfacial tension (IFT) between the phases and θ ij is the contact angle that the main terminal meniscus (MTM) makes with the solid surface. For pores with angular cross section, in addition to MTM, an arch meniscus (AM) forms at each corner. Assuming a pore with straight walls, the capillary pressure across the AM can be obtained from the Young-Laplace equation Pcthij = Pi − Pj =

γij rij

,

(2)

where rij is the radius of the arch forming at the corners between phases i and j. The threshold capillary pressure for displacement in pores with angular cross sections can be formulated using the MS-P method, which is based on the minimization of the Helmholtz free energy when a wetting phase is displaced by a non-wetting phase. The details of this method are reported elsewhere.24 For the oil-water drainage inside a pore with an angular cross section, the threshold capillary pressure is derived as Pcthow =

γow [Low,t + Los,t cos(θow )] , Ao,t

(3)

where Lij,t is the total length along which phases i and j are in contact with each other and Ao,t is the total area of the cross section that oil occupies (here, subscripts “w,” “o,” and “s” denote water, oil, and solid, respectively). For a square cross section with a perimeter of Lt (subscripts “t” and “c” denote total and corner, respectively), we have Low,t = 4Low,c , Los,t = Lt − 8b, and Ao,t = At − 4Aw,c , where Low,c , b, and Aw,c are the length of contact line between oil and water, meniscusapex distance of the interface, and the area occupied by water at a corner, respectively. With simple trigonometric algebra it can be shown that π − (θow + α) , (4) Low,c = 2row 2 b = row

cos(θow + α) , sin(α)

(5)


194703-3


J. Chem. Phys. 141, 194703 (2014)

π cos(θow + α) Aw,c = cos(θow ) − − (θow + α) , sin(α) 2 (6) where α is the corner half angle. Equation (3) can be simplified to Eq. (1) for the case of circular cross section. The threshold capillary pressure can also be formulated based on the same principle when the angular pore is mixedwet with water in the corners and oil in the center and we displace the oil with water in a piston-like displacement. This process is called waterflooding and can produce different pore fluid occupancies depending on the wettability of the pore. In water-wet pores, it leads to a pore fully saturated with water, whereas in oil-wet pores, oil may leave a wetting layer between water in the corners and in the center depending on the oil-wetness of the pore and stability of the layer. The stability of the oil layer can be examined based on the threshold capillary pressure using the following approach. If the threshold capillary pressure with an oil layer is smaller than the one without the oil layer, then the layer is stable as long as the capillary pressure remains between those two thresholds and it starts to collapse as the capillary pressure becomes equal to or greater than the second threshold. However, if the threshold capillary pressure for the latter case is in fact smaller than that for the former case, then the formation of an oil layer is not thermodynamically feasible and this layer will not form at the beginning of waterflooding. We refer the interested readers to Refs. 35 and 36 where the threshold capillary pressure formulations for every possible scenario in waterflooding are presented with details. Here, we only list the calculation results with the threshold capillary pressure equations for waterflooding. 2 row

B. Molecular dynamics simulation

In this work, the coarse grained (CG) models suggested by Marrink et al.22, 23 were employed to represent water and oil molecules and wall particles. In these CG models, on average, every four atoms (excluding hydrogens) are represented by one coarse grained particle. Therefore, every four water molecules form one bead with a mass of 72 amu and dodecane molecules, representing oil in our simulations, consist of 3 beads, with a mass of 56.7 amu, attached to each other with the harmonic bond potential given by 1 K (r − σ )2 , 2 bond ij

(7)

where Kbond is the bond force constant, which in our simulations is equal to 1250 kJ/(mol nm2 ). The angle between the coarse grains in the oil molecules was kept flexible with a force constant of 25 kJ/mol using the following angle potential:

Uang (θ ) =

1 K (cos(θ ) − 1)2 . 2 ang

Water Oil

Water

Oil

Wall1 (SWW)

Wall2 (WOW)

Wall3 (SOW)

5.0 2.0

2.0 3.4

5.0 2.0

2.0 3.1

2.0 4.5

The intermolecular forces are applied using shifted Lennard Jones potentials as shown below: ULJ (r) − ULJ (rc ) + FLJ (rc )(r − rc ) r < rc ULJ (r) = 0 r > rc (9) with 6

σij 12 σij − ULJ (r) = 4εij (10) r r and FLJ (r) = −

∂ULJ (r) . ∂r

(11)

The sigma of all particles was equal to 0.47 nm and the cut-off radius was set to 1.5 nm. The LJ potentials for the fluids and pore wall particles are reported in Table I. Three wettabilities were considered for the walls, which are labeled as strongly water-wet (SWW), weakly oil-wet (WOW), and strongly oil-wet (SOW). The piston in our MD simulation systems is a flat wall consisting of hard sphere particles allowed to move in only one direction (i.e., Z-direction), which is controlled by the external pressure exerted on them. 2. Nanopore structures

1. Force field parameters

Ubond (rij ) =

TABLE I. LJ potentials (in kJ/mol) of the particles used in the MD simulations.

(8)

We configured certain nanopore structures for MD simulations in this study. The wall particles in nanopores have nonordered configurations with a density of 0.846 σ −3 . Three cross-sectional shapes were considered for the nanopores: circle, square, and right triangle (with acute angles of 30◦ and 60◦ ), as illustrated in Figure 1. The pores have a height of 10 nm and similar cross section areas varying between 180 and 182 nm2 . The thickness of the pore walls at any position is always equal to or higher than 1.5 nm, so that fluid particles that are very close to the walls cannot feel the empty space behind them. Figure 1 shows circular cross sections with smooth and rough walls. To randomize the roughness and avoid creating sharp and unrealistic indentations, the wall particles were divided into inner and outer sections. We only allow the inner particles to move freely while randomly assigning different interaction potentials between the inner and outer particles. The rough wall was only considered for the pore with circular cross section in drainage simulations (which is explained in later sections), and the roughness was calculated from the equation below25 as described in Ref. 4, N 2 i=1 (di − di ) , (12) R= N −1


194703-4


J. Chem. Phys. 141, 194703 (2014)

FIG. 1. Pore structures used in the MD simulations: (a) Circle with smooth wall; (b) circle with rough wall; (c) square; and (d) right triangle.

where di is the minimum distance between the longitudinal axis of the pore and the wall particles in the ith bin, di is the average distance, and N is the total number of bins. The values of roughness for the smooth and rough walls were found to be 0.06 and 0.7, respectively. In the MD simulation studies mentioned earlier, the wall particles are modeled as either frozen or restrained to their equilibrium positions with a spring force or a finitely extensible nonlinear elastic (FENE) potential.4–16 In this work, a few preliminary MD simulations were run to determine if thermal vibration of pore walls can significantly affect the simulation outcomes. For the drainage tests in the circular pore, MD simulations were performed with frozen and tethered wall particles using the harmonic potentials with force constants of 1000, 3000, and 5000 kJ/(mol nm2 ). Except for the simulations with the force constant of 1000 kJ/(mol nm2 ) where there was a considerable amount of fluids escaping through the walls, the trajectories of the simulations were similar. The invasion of oil into the pore was slightly faster for frozen particles than restrained ones. A similar conclusion was reached in Ref. 4. That is, the thermal agitation of wall particles did

not significantly impact the dynamics of dodecane molecules invading into capillaries during spontaneous imbibition simulations. In our simulations, we have used frozen walls, since the diffusion of fluid particles through the wall is the lowest, and therefore we were able to use a minimum wall thickness of 1.5 nm. 3. MD parameters

Generally, NVE (constant number of particles, volume, and energy) ensembles can produce the true hydrodynamics of the system of study; however, for nonequilibrium molecular dynamics (NEMD) simulations performed by applying an external pressure on the piston, the temperature would rise as a result of the external energy that is being pumped into the system. Therefore, a thermostat is needed to dissipate this heat. The choice of thermostat is important for these NEMD simulations since we have center of mass movements. The thermostats that use the global kinetic energy to control temperature are susceptible to the flying ice cube effect.26 Among the thermostats that control temperature locally, dissipative


194703-5


J. Chem. Phys. 141, 194703 (2014)

particle dynamics (DPD) approach is shown to preserve the momentum of the system and is preferred over Langevin thermostat when the hydrodynamic properties of the system are important.27 In DPD thermostat, in addition to the Van der Waals (VDW) force (Fi ), a dissipative force (FiD ), and a random force (FiR ) are applied to each particle,28 p˙i = Fi + FiD + FiR .

(13)

The dissipative force is designed to act as a drag force and depends on the relative velocity ( vij ) between the particles. It is calculated from this equation, FijD = −ξ w 2 (rij )(rij · vij )ˆrij ,

(14)

where ξ is the friction factor, which in our simulations was set to 0.5 (gr/mol ps), and w(r) is the weight function defined as w(r) = 1 − r/rc for r < rc and w(r) = 0 for r ≥ rc . The cutoff radius (rc ) is set to 1.5 nm. The random forces are given by FijR = ς w(rij )αt −1/2 rˆij ,

(15)

where ς = (2ξ kT)0.5 and is designed to represent a stochastic heat bath, and α is a Gaussian random number with zero mean and unit variance. An alternative approach to stochastic thermostats is to only thermostat the wall particles. This method has the advantage of not perturbing the dynamics of fluid particles since there is no thermostat term in their Hamiltonian equations.29 However, we did not consider this technique in our simulations since we are looking at the final equilibrium state of fluids. Instead, we used the DPD thermostat so that we could freeze the wall particles in order to speed up the simulations and reduce the diffusion of fluid particles out of the pore. Also the stochastic forces in this thermostat helped prevent CG water from freezing inside the pore. The time step was chosen to be 0.005 τ , where τ = (σ 2 m/kB T)0.5 . Periodic boundary conditions (pbc) were used in all directions, however, we increased the length of the simulation box in the Z-direction so that the minimum distance between the displaced fluid surface and the piston (or its pbc image) is at least 4 nm in order for the displaced fluid to be in equilibrium with its vapor phase in the vacuum region between its surface and the piston. The MD simulations were run using DL_POLY_4.0.2 simulation package.30 The source code was modified to implement the DPD thermostat and corrected for the piston displacement in the External Field section. The MD snapshot pictures were produced with VMD software.31 4. Interfacial tension and contact angle

The IFT and the contact angle used in the MS-P method to determine threshold capillary pressures, were obtained from MD simulations as explained in this section. The IFT between water and oil was calculated by placing slabs of water and oil with a size of 10 × 10 × 5 nm3 in contact with each other at a temperature of 302 K and different pressures. At each pressure, an NPT (constant number of particle, pressure, and temperature) ensemble was run for 4 ns by employing the

Nose-Hoover thermostat and barostat. IFT values were calculated during the last 3 ns using the pressure tensor components in the following equation:32

Lz 1 Pzz − (Pxx + Pyy ) , γ = (16) 2 2 where γ is the interfacial tension, Lz is the simulation box size in the Z-direction, which is normal to the oil and water surface, and Pii is the ensemble average of normal pressure in the i direction. To check if using the DPD thermostat would change the IFT values, NVT (constant number of particle, volume, and temperature) ensembles were performed for 4 ns with the DPD thermostat following the last configuration of the NPT simulations. The IFT values from both simulations were very similar. The values from NPT simulations were 49.8, 50.8, and 50.2 mN/m for pressure values of 100, 125, and 150 bar, respectively. These values are in the range of experimentally measured IFT counterparts for water and alkanes.33 Although the IFT values should increase with pressure,34 the increase is negligible and is within the simulation margin of error. The equilibrium contact angle was determined in the static conditions by inserting a cubic box of the wetting fluid (either water or oil depending on the wettability of the wall) with a side of 5 nm inside a larger box of the non-wetting phase with a size of 20 nm × 20 nm × 10 nm, placed over a slab of wall particles. The NVT simulations were run for 4 ns using DPD thermostat at a temperature of 302 K. The contact angle was obtained from the last snapshot of MD simulations using ImageJ software.37 Figure 2 displays snapshots of the initial and final configurations of the MD simulation performed to measure the equilibrium contact angle on the strongly water-wet wall. The size dependency of IFT and contact angle was examined by running similar simulations at a larger scale. For IFT calculation, a simulation box of 20 × 20 × 40 nm3 was prepared by placing two cubic simulation boxes of water and dodecane molecules adjacent to each other in the Z-direction. The IFT obtained from an NPT simulation at 125 bar was 50.5 mN/m, which is in the range of IFT values obtained from the small-scale simulation. To determine the contact angle in a water-wet pore, a cubic box of (7.5 nm)3 full of water was immersed into dodecane molecules confined inside a rectangular box of 36 × 36 × 16 nm3 that was placed on a slab of wall particles with a thickness of 2 nm. These dimensions are approximately 3 times larger than those used in the smallscale simulation. An NPT simulation was run for 9 ns and the calculated contact angle was in the range of 15◦ –20◦ , in good agreement with the small-scale simulation. Similar to Stukan et al.,4 the dynamic (i.e., advancing and receding) contact angles in the smooth circular pore were calculated from the shape of the meniscus. The pore space was divided into 10 concentric cylindrical shells with equal cross sectional area and the fluid level was determined at each shell. The heights of these levels were then plotted versus the distance of the shell from the central axis of the pore and a circular curve was fitted through the data. The slope of the tangent line to the curve at a point located at a distance of 0.5 nm from the wall, represents the dynamic contact angle. For the


194703-6


J. Chem. Phys. 141, 194703 (2014)

(a)

(b)

FIG. 3. MD simulation of oil spontaneous imbibition in a nanopore: (a) Meniscus position versus time; and (b) dynamic contact angle versus time.

FIG. 2. Snapshots of the MD simulation performed to measure the equilibrium contact angle on the strongly water-wet wall: Oil, water, and wall particles are colored as red, pink, and green, respectively; (a) Initial configuration and (b) final configuration.

spontaneous imbibition simulations (Figure 3), the dynamic contact angles were then fitted to the following equation in order to determine the equilibrium contact angle:4 cos θ =

l(t) cos θe , l(t) + P2 /P3

(17)

where θ and θ e are the dynamic and equilibrium contact angles, respectively, and l (t) is the distance from the meniscus to the nanopore entry. P2 and P3 are fitting parameters in the modified Lucas-Washburn equation which are further elaborated in Ref. 4. III. RESULTS AND DISCUSSION

To test the validity of our simulation setups, we initially simulated the spontaneous imbibition of oil molecules in oil-wet cylindrical nanotubes with a radius of 7.6 nm and a length of 47 nm. The wall had a smooth surface and two levels of wettabilities were considered: εo−pw = 1.4 kT

= 3.49 kJ/mol (strongly oil-wet) and εo−pw = 1.1 kT = 2.74 kJ/mol (weakly oil-wet). Similar to previous studies,4, 6–8 we found that the dynamic of oil rise in the nanotube follows the modified Lucas-Washburn equation. This indicates that the hydrodynamics of fluid flow in our simulations are modeled correctly. Figure 3 displays the change in the meniscus level and the dynamic contact angle versus time. The equilibrium contact angle was determined from two methods: (1) by fitting Eq. (17) to the dynamic contact angles calculated during the simulation (the dashed lines in Figure 3(b)), and (2) by placing a cubic box of oil molecules on a slab of wall particles, as explained in Sec. II. The contact angles calculated from the first method were 15.2◦ and 71.5◦ and those measured using the second method were in the range of 15◦ -20◦ and 60◦ –70◦ for strongly and weakly oil-wet nanotubes, respectively. The values show reasonable agreement.

A. Drainage

For the oil-to-water drainage simulations, we initially packed water molecules inside the pore and placed two slabs of water, with 2 nm thickness each, at the entrance and exit of the pore. We then placed dodecane molecules, with a volume 2.2 times greater than the pore volume, between the lower water slab and the piston. An initial 200 ps of equilibrium simulation was run before applying any forces on the piston. The initial configuration for the circular tube is depicted in Figure 4.


194703-7


J. Chem. Phys. 141, 194703 (2014)

FIG. 5. Calculated external pressure from an NVE simulation using a cylindrical frozen wall.

FIG. 4. Initial configuration of drainage simulation for a circular pore.

The first simulation was performed to verify if the external pressure defined in DL_POLY is the actual pressure that is exerted on the system. To this end, we considered the NVE ensemble (i.e., no thermostat) so that the change in total energy of the system is equal to the work that is being done, as there is no heat transfer (U = W + Q, Q = 0). This work is equal to the product of pressure by the change in volume so that W = Pext V = Pext Az. Therefore, we can determine the external pressure that is being applied on the system at any time interval from the following relation: Pext = U/Az,

(18)

where A is the area of the box in the X-Y plane and z is the position of the piston in the Z-direction. Different NVE simulations were performed for both circular and square cross sections while the external pressure was set to 80, 100, and 120 bar. However, the external pressure that was calculated from Eq. (18) for each simulation was greater than the preset pressure and was dependent on the domain decomposition used in the system. We modified the source code to correct for this error. We then tested the correct code for preset external pressures of 100, 120, and 140 bar for both frozen and tethered wall particles and found that the calculated external pressures matched the preset pressures within 5% in all simulations. The calculated external pressure during simulations with a cylindrical frozen pore and the preset external pressure of 120 bar are shown in Figure 5. The figure shows the build-up in the external pressure as the piston pushes oil until it reaches the preset value. The threshold capillary pressure can be calculated for the pores with smooth walls in the three cross sections shown in Figure 1 using the MS-P method (i.e., Eqs. (1)–(6)) for which IFT and contact angle are input parameters. As was mentioned in Sec. II, the IFT values obtained from MD simulations were

in the range of 50–51 mN/m and the equilibrium contact angles for the strongly water-wet wall were in the range of 15◦ – 20◦ . The equilibrium contact angle was used for the MS-P calculations since, as it will be discussed later, advancing and receding contact angles in our systems of study were similar and equal to the equilibrium contact angle. Using the IFT value of 50.5 mN/m and the contact angle of 17.5◦ , the values of Pcthow were calculated for the pores with circular, square, and triangular cross sections. The results are presented in Figure 6. To check the validity of the MS-P method in nanopores, we performed drainage simulations with variable external pressures. In the NEMD simulations, oil starts to invade the pore at a pressure higher than the value predicted by the MS-P method. After oil invades the pore, we decreased the external pressure by 5 bar every 1 ns until the oil begins to retract from the pore. The pressure at which the oil can no longer remain inside the pore is the threshold capillary pressure. At the end of each simulation, snapshots of the system taken every 200 ps were analyzed and the number of oil and water particles inside the pore was determined. The profile of oil and water particles inside the pore versus simulation time for the circular cross section is shown in Figure 7(a). As illustrated in this figure, oil particles inside the pore are replaced by water particles as the external pressure changes from 135 to 130 bar (the fact that the number of water and oil particles matches each other is purely coincidental). Considering that the

FIG. 6. Threshold capillary pressures obtained from the MD simulations and the MS-P method for drainage in the nanopores with different cross-sectional shapes.


194703-8


J. Chem. Phys. 141, 194703 (2014)

(a)

(b)

(c)

(d)

FIG. 7. Number of particles inside the pore versus simulation time: (a) Number of oil and water particles for the circular cross section; (b) number of oil particles for the three cross-sectional shapes with smooth wall; (c) number of oil particles inside the circular pore with smooth wall when the external pressure drops from 150 to 0 bar after 2 ns; and (d) number of oil particles inside the pore with square cross section.

water pressure is equal to its vapor pressure at 302 K, which is much smaller than 1 bar (∼0.04 bar), it is safe to assume that the pressure difference between oil and water is equal to the external pressure on the piston. Hence, the threshold capillary pressure for the circular cross section is in the range of 135–130 bar. Figure 7(b) displays the number of oil particles inside the pore versus simulation time, for the three cross sections. From these profiles, we can determine the ranges of Pcthow for each pore. The Pcthow obtained from the MD simulations are compared with the MS-P counterparts in Figure 6. The results show an encouraging agreement. For the circular pore, MD simulations predict a slightly higher threshold capillary pressure than the MS-P method. This could be due to the fact that a thin film of water remains on the wall surface during drainage, which reduces the pore cross-sectional area available to oil. This water film is less significant for the pores with angular cross sections since almost all water molecules retract to the corners. One question to consider here is how fast the oil particles inside the pore would feel the change in the external pressure that is being applied on the piston. If the latency time is more than 1 ns, then the threshold capillary pressure would be higher than the value predicted by the simulations. To answer this question, the delay in the response of the meniscus

to a change in external pressure was examined by performing a simulation where the external pressure dropped from 150 to 0 bar after 2 ns. The results shown in Figure 7(c) indicate that the response of oil particles to the external pressure is almost instantaneous and thus 1 ns is a sufficient amount of time for oil to react to the new pressure exerted on the piston. We also investigated the advancing and receding contact angles. Receding contact angle is measured when the nonwetting phase is displacing the wetting phase during drainage and is usually smaller than the advancing contact angle measured during imbibition. This is called contact angle hysteresis and is attributed to wall roughness, change of wall wettability, etc. Contact angle hysteresis results in different threshold capillary pressures depending on whether the nonwetting phase is invading or receding. To check the variation of capillary pressure, we ran drainage MD simulations for the square cross section at external pressures of 135 and 145 bar. The number of oil particles inside the pore during these simulations is depicted in Figure 7(d). The results suggest that Pcthow is almost the same whether the oil is invading or retrieving, meaning that the advancing and receding contact angle must be the same. Similar results were obtained for other cross sections as well as the circular pore with the rough wall. This may suggest that the wall roughness does not have a noticeable


194703-9


(a)

J. Chem. Phys. 141, 194703 (2014)

the smooth and rough surfaces was similar (i.e., 17%), although the remaining water distributions were different. The main difference between these pores was that a water film uniformly covered the surface of the smooth pore, while for the rough wall, water stayed in the grooves, as shown in Figures 9(a) and 9(b). The configurations of oil and water after drainage in pores with different cross sections are illustrated in Figures 9(a)–9(d). As expected, the amount of water that resides in the corners increases as the corner angle decreases. The residual water saturation at the end of drainage is 17%, 20%, and 36% for pores with circular, square, and triangular cross sections, respectively. The fact that the threshold capillary pressures obtained from MD simulations are consistent with the values calculated from Eqs. (1)–(6), implies that the interfacial tension between oil and water remained almost constant inside the pore. In nanoscale pores, the interfacial tension between fluids may change inside the pore due to wall effects. However, in our simulations, we used the short range interactions (limited to 1.5 nm), and the density of the wall particles and their interaction potentials were similar to those of fluids. As future work, we plan to add polar and charged particles to the simulation system in order to include long range electrostatic interactions and examine their effects on the results.

(b)

B. Waterflooding FIG. 8. Drainage simulation in the smooth circular pore: (a) Advancing and receding meniscus at about 7 nm from the pore entrance; and (b) dynamic contact angles versus time.

impact on the capillary pressure and contact angle hysteresis at this scale. For a circular cross section, the receding and advancing contact angles can be calculated from the shape of the meniscus using the method that was described in Sec. II. Due to the symmetry of oil-water interface in the circular pore, only half of the advancing and receding meniscus at a same position (the meniscus level is about 7 nm from the pore entrance), are shown in Figure 8(a). The dynamic contact angle was estimated from the circular curve at a distance of 0.5 nm from the wall since the wall is covered with one layer of water molecules. Figure 8(b) shows the receding and advancing contact angles during the simulation. The contact angle starts to decrease as the curvature of the oil-water interface increases until it reaches the value of the receding contact angle. Subsequently the meniscus starts to move with a constant curvature as the contact angle remains fixed. Figure 8(b) also shows that the advancing and receding contact angles have similar values (i.e., 18◦ –20◦ ). This is consistent with the equilibrium contact angles obtained from the MD simulations and therefore we used equilibrium value in the MS-P calculations. These results suggest that in our pore systems the effect of wall roughness on Pcthow was not significant as the radius of the cross section did not significantly vary inside the pore. Additionally, the remaining water saturation in the pores with

After injecting more than 2 PV of oil into each pore, we used a snapshot of the pore structure that included 1.5 nm of the medium below the pore entrance and above the pore exit, and placed a slab of oil with a thickness of 1 nm on both sides. Similar to drainage, 2.2 PV of water was placed between the lower oil slab and the piston. The wettability of wall particles were changed based on the configuration of oil and water inside the pore at the end of drainage. As previously mentioned, the wall surface in the corners of angular cross sections retained their wettability since they did not come in contact with the non-wetting phase. For the circular cross section, the wettability of all wall particles was altered to strongly oil-wet. For the angular cross sections, the wall particles that came in contact with oil experienced wettability alteration to strongly oil-wet. For the square and the triangular cross sections, the percentages of oil-wet wall particles were 51% and 53%, respectively. Figures 9(c) and 9(d) illustrate the change of wall wettability for the angular cross sections (the cyan particles are oil-wet). For the pore with square cross section, in addition to strongly oil-wet wettability, we also considered weakly oil-wet wall particles with smaller εoil−wall value as reported in Table I. For waterflooding simulations, the threshold capillary pressure was first determined with the same method used for drainage simulations; i.e., starting with high external pressures then decreasing pressure incrementally by 10 bar every 1 ns until the invading water starts to retreat. The threshold capillary pressure was also estimated from the MS-P method using an IFT value of 50.5 mN/m and an equilibrium contact angle of 150◦ and 170◦ for the weakly and strongly oil-wet cases, respectively (as shown previously in Figure 8,


194703-10


J. Chem. Phys. 141, 194703 (2014)

FIG. 9. Final configurations of oil (red particles) and water (pink particles) in pores with different cross sections and wall roughness (water-wet and oil-wet wall particles are colored green and cyan, respectively): (a) Circle with smooth wall; (b) circle with rough wall; (c) square; and (d) triangle.

advancing and receding contact angles were equal to the equilibrium contact angle). While we acknowledge that weakly oil-wet may refer to contact angle of less than 135◦ , the notions of weakly and strongly oil-wet are applied merely to distinguish the levels of wettability alterations used in our MD simulations. In Figure 10, the results of the MS-P calculations are compared with the results of MD simulations. For the circular cross section, the threshold capillary pressure is smaller than the MS-P calculation by about 10%. The lower threshold capillary pressure might be attributed to the water film that initially covered the wall surface. The water film makes it easier for the water to invade since the invading water particles are in contact with water instead of wall surface. The thickness of this water film is less than a few nanometers, therefore its effect on the threshold capillary pressure would only be important for the nanoscale pore structures. For the angular cross sections, there is a very good agreement between the MD simulations and the MS-P results, which suggests that MS-P method is still valid for waterflooding in nanoscale confinements consisting of Lennard-Jones particles.

The waterflooding simulations were performed at the threshold capillary pressure so that flow would be capillary dominated. During the simulations in the angular pores, water invaded from the center of the pores and therefore, for

FIG. 10. Threshold capillary pressures obtained from the MD simulations and the MS-P method for waterflooding in nanopores with different crosssectional shapes.


194703-11


J. Chem. Phys. 141, 194703 (2014)

continued. The final configuration and the oil layer disruption are displayed in Figure 11 for the triangular cross section (with strongly oil-wet wall particles). The instability of the oil layer was observed in both triangular and square cross sections with different oil-wet walls. From the MS-P method, the threshold capillary pressure for the scenarios with oil layer formation in the square and triangular cross sections with strongly oil-wet walls, were 139.8 and 117.8 bar, respectively, and in the square with weakly oil-wet walls, it was 126.6 bar. The fact that the threshold capillary pressure for the pore having an oil layer is higher (in absolute value) than that of the pore without any oil layer means that the oil layer is not thermodynamically stable, which is in line with our MD simulations. It should be noted that as the corner half angle decreased, the stability of the oil layer increased and remained stable for a longer time. After water displaced oil inside the nanopores, thin films of oil remained on the oil-wet parts of the walls (as shown in Figure 11(c)). IV. CONCLUSIONS

Molecular dynamics simulations on large parallel computers were used to investigate some of the key aspects of multiphase fluid/fluid displacements in nanoscale pores with different cross sectional shapes (i.e., circle, square, and right triangle). Each nanopore structure included an impermeable piston assembly to enable forced fluid/fluid displacements in the pore. The DL_POLY source code was corrected so that the pressure felt by the pore is equal to the external pressure assigned to the impermeable piston. The drainage simulations were performed by injection of water into oil in water saturated, water-wet pores. The threshold capillary pressures obtained from the MD simulations were in good agreement with those calculated using the MS-P method. The order of Pcthow and the remaining water saturation in pores with different cross-sectional shape, but same cross-sectional area, was as follows: circle < square < right triangle. The nanopores with angular cross sections retained the wetting phase (water) in the corners at the end of drainage. Following oil drainage, we allowed wall particles contacted by oil to experience wettability alteration. We investigated the impact of wettability alteration on the subsequent waterflooding and the resulting pore fluid configurations. For the advancing and receding water, the Pcthow and the contact angles were the same, within the precision of our MD simulations. During waterflooding simulations in oil-wet pores with angular cross sections, oil layers were formed over water surface in the corners, although, they were not stable and collapsed as water continued to invade the pore. The instability of oil layers was also predicted from the MS-P method for the conditions of MD simulations. The oil-wet walls however, retained a film of oil particles. FIG. 11. Final configuration and oil layer disruption in a pore with triangular cross section: (a) Oil layer is formed; (b) oil layer collapsed; and (c) final configuration.

ACKNOWLEDGMENTS

a short period of time there was an oil layer between the water in the center and the corner. However, this oil layer was not stable and eventually collapsed as the simulation

We gratefully acknowledge financial support of Hess Corporation and the school of Energy Resources at the University of Wyoming. We thank Dr. Will Welch and Dr. Jan Kubelka for many insightful discussions on this work. We


194703-12


also thank Dr. Arsalan Zolfaghari for helping with the MSP method calculations. 1 H.

Darabi, A. Ettehad, F. Javadpour, and K. Sepehrnoori, J. Fluid Mech. 710, 641 (2012). 2 E. Fathi and I. Y. Akkutlu, SPE J. 18, 27 (2013). 3 V. P. Sokhan, D. Nicholson, and N. Quirke, J. Chem. Phys. 115, 3878 (2001). 4 M. R. Stukan, P. Ligneul, J. P. Crawshaw, and E. S. Boek, Langmuir 26, 13342 (2010). 5 W. Stroberg, S. Keten, and W. K. Liu, Langmuir 28, 14488 (2012). 6 G. Martic, F. Gentner, D. Seveno, D. Coulon, J. De Coninck, and T. D. Blake, Langmuir 18, 7971 (2002). 7 G. Martic, F. Gentner, D. Seveno, J. De Coninck, and T. D. Blake, J. Colloid Interface Sci. 270, 171 (2004). 8 G. Martic, T. D. Blake, and J. De Coninck, Langmuir 21, 11201 (2005). 9 D. I. Dimitrov, A. Milchev, and K. Binder, Phys. Rev. Lett. 99, 054501 (2007). 10 D. I. Dimitrov, A. Milchev, and K. Binder, Langmuir 24, 1232 (2008). 11 S. Gruener and P. Huber, Phys. Rev. Lett. 103, 174501 (2009). 12 S. Supple and N. Quirke, J. Chem. Phys. 121, 8571 (2004). 13 S. Supple and N. Quirke, J. Chem. Phys. 122, 104706 (2005). 14 C. Chen, C. Gao, L. Zhuang, X. Li, P. Wu, J. Dong, and J. Lu, Langmuir 26, 9533 (2010). 15 C. Chen, L. Zhuang, X. Li, J. Dong, and J. Lu, Langmuir 28, 1330 (2012). 16 L. Joly, J. Chem. Phys. 135, 214705 (2011). 17 A. R. Kovscek, H. Wong, and C. J. Radke, AIChE J. 39, 1072 (1993). 18 M.-H. Hui and M. J. Blunt, J. Phys. Chem. B 104, 3833 (2000). 19 M. Piri and M. J. Blunt, Phys. Rev. E 70, 061603 (2004). 20 M. Piri and M. J. Blunt, Phys. Rev. E 71, 026301 (2005). 21 S. Joseph and N. R. Aluru, Nano Lett. 8, 452 (2008).

J. Chem. Phys. 141, 194703 (2014) 22 S.

J. Marrink, A. H. de Vries, and A. E. Mark, J. Phys. Chem. B 108, 750 (2004). 23 S. J. Marrink, H. J. Risselada, S. Yefimov, D. P. Tieleman, and A. H. de Vries, J. Phys. Chem. B 111, 7812 (2007). 24 Z. T. Karpyn and M. Piri, Phys. Rev. E 76, 016315 (2007). 25 A. L. Sumner, E. J. Menke, Y. Dubowski, J. T. Newberg, R. M. Penner, J. C. Hemminger, L. M. Wingen, T. Brauers, and B. J. Finlayson-Pitts, Phys. Chem. Chem. Phys. 6, 604 (2004). 26 J. E. Basconi and M. R. Shirts, J. Chem. Theory Comput. 9, 2887 (2013). 27 C. Pastorino, T. Kreer, M. Müller, and K. Binder, Phys. Rev. E 76, 026706 (2007). 28 T. Soddemann, B. Dünweg, and K. Kremer, Phys. Rev. E 68, 046702 (2003). 29 H. Frentrup, C. Avendao, M. Horsch, A. Salih, and E. A. Müller, Mol. Simul. 38, 540 (2011). 30 I. T. Todorov, W. Smith, K. Trachenko, and M. T. Dove, J. Mater. Chem. 16, 1911 (2006). 31 W. Humphrey, A. Dalke, and K. Schulten, J. Mol. Graph. 14, 33 (1996). 32 G. J. Gloor, G. Jackson, F. J. Blas, and E. de Miguel, J. Chem. Phys. 123, 134703 (2005). 33 S. Zeppieri, J. Rodrguez, and A. L. Lpez de Ramos, J. Chem. Eng. Data 46, 1086 (2001). 34 M. Abdelrahim, “Measurement of interfacial tension in hydrocarbon/water/dispersant systems at deepwater conditions,” M.S. thesis (Louisiana State University, 2008). 35 A. Z. Shahrak, “Pore-scale network modeling of two- and three-phase flow based on thermodynamically consistent threshold capillary pressures,” Ph.D. dissertation (University of Wyoming, 2014). 36 A. Z. Shahrak and M. Piri, “Pore-scale network modeling of two- and three-phase flow based on thermodynamically consistent threshold capillary pressures” (unpublished). 37 C. A. Schneider, W. S. Rasband, and K. W. Eliceiri, Nat. Methods 9, 671 (2012).


Molecular Dynamics Simulation of the Effect of Angle Variation on Water Permeability through Hourglass-Shaped Nanopores.

Molecular simulation of water vapor outgassing from silica nanopores.

Sensing of protein molecules through nanopores: a molecular dynamics study.

Wetting of graphene oxide: a molecular dynamics study.

Presence of electrolyte promotes wetting and hydrophobic gating in nanopores with residual surface charges.

Ar matrices: molecular dynamics and luminescence study.

Molecular dynamics analysis on wetting and interfacial properties of water-alcohol mixture droplets on a solid surface.

Study on the characteristics of gas molecular mean free pathin nanopores by molecular dynamics simulations.

Molecular dynamics study on the wettability of a hydrophobic surface textured with nanoscale pillars.

Influence of water content in mixed solvent on surface morphology, wettability, and photoconductivity of ZnO thin films.

Effects of surface wettability and liquid viscosity on the dynamic wetting of individual drops.

Energetics and Predissociation Dynamics of Small Water, HCl, and Mixed HCl-Water Clusters.

Coupling all-atom molecular dynamics simulations of ions in water with Brownian dynamics.

Water Wettability.

Origin of the Shape of Current-Voltage Curve through Nanopores: A Molecular Dynamics Study.

Molecular dynamics study of wetting behavior of grafted thermo-responsive PNIPAAm brushes.

Water electrolyte transport through corrugated carbon nanopores.

Wetting transitions of water on graphite and graphene.

Molecular dynamics simulation of interlayer water embedded in phospholipid bilayer.

Stability and dynamics of membrane-spanning DNA nanopores.

Deswelling of ultrathin molecular layer-by-layer polyamide water desalination membranes.

Molecular dynamics simulations of the nucleation of water: determining the sticking probability and formation energy of a cluster.

A general and facile chemical avenue for the controlled and extreme regulation of water wettability in air and oil wettability under water.

Electric double-layer structure in primitive model electrolytes: comparing molecular dynamics with local-density approximations.