Molecular dynamic study of capillary forces on clay particles

Applied Clay Science 88–89 (2014) 170–177

Contents lists available at ScienceDirect

Applied Clay Science journal homepage: www.elsevier.com/locate/clay

Research paper

Molecular dynamic study of capillary forces on clay particles Priyanthi M. Amarasinghe a, A. Anandarajah a,⁎, Pijush Ghosh b a b

Department of Civil Engineering, Johns Hopkins University, MD, USA Department of Applied Mechanics, Indian Institute of Technology, Madras, India

a r t i c l e

i n f o

Article history: Received 3 February 2012 Received in revised form 16 August 2013 Accepted 13 December 2013 Available online 15 January 2014 Keywords: Capillary Contact angle Surface tension Molecular dynamics Pyrophyllite

a b s t r a c t Capillary forces between discrete colloid particles are of interest in many areas of science and engineering. This study used molecular dynamics (MD) to examine capillary phenomena that develop in a clay–water–air system, which requires an understanding of capillary forces, contact angles, and meniscus curvatures. The capillary formed between two parallel pyrophyllite clay particles and pure water was studied using MD for various system parameters which included particle separation, particle thickness, and strength of the clay–water van der Waals force. For a given set of potential energy parameters, the contact angle was independent of particle separation (≈ 550), the radius of curvature depended on particle separation, and particle thickness had a negligible influence on contact angles and meniscus radii. The van der Waals energy parameters, however, had a profound influence on contact angles and meniscus radii. In all cases, capillary forces on clay particles calculated by MD agreed with values calculated using the classical Young–Laplace equation. When the separation distance is 40 Å, the force was calculated to be 9.87 × 10−10 N by Young–Laplace equation and 8.37 × 10−10 N by MD. © 2013 Published by Elsevier B.V.

1. Introduction Capillary effects are of interest in many fields such as soil science, material science, colloidal science and engineering. For instance, when a dry soil becomes partially saturated, capillary forces develop between particles, resulting in shrinkage in some soils such as clays. Numerous past studies have employed molecular dynamics to further understand the phenomenon of capillarity. Previous studies, which examined capillary flow in narrow pores, used structureless walls or arranged atoms on a lattice to simulate a solid wall (Bitsanis et al., 1987; Heinbuch and Fischer, 1989; Koplik et al., 1989; Martic et al., 2002; Sokołowski, 1991; Thompson and Robbins, 1989; Trozzi and Ciccotti, 1984). De Ruijter et al. (1999) used molecular dynamics to study spreading of liquid droplet on a flat solid surface. The solid surface was modeled using atoms placed on a face-centered cubic lattice. Molecular dynamic technique is proven to be an effective method in studying clay structures and their behaviors. For instance, Nakano et al. (2003) studied adsorption of cesium on clay surfaces. Ichikawa et al. (2004) studied the diffusion of HTO (tritium water) in clay structure using molecular dynamics. The pair-wise interaction between atoms was simulated using the Lennard–Jones (LJ) 12–6 potentials. The contact angle significantly changed as the solid–liquid LJ potential parameters were varied. Martic et al. (2002) showed that contact angle is dependent on the system wetting velocity. There have been numerous attempts to use experiments to understand the surface tension of clay mineral–surfactant interfaces (Cipriano et al., 2005; Docoslis et al., 2000; Lee et al., 2000; ⁎ Corresponding author at: Department of Civil Engineering, Johns Hopkins University, MD 21218, USA. Fax: +1 410 516 7473. E-mail address: [email protected] (A. Anandarajah). 0169-1317/$ – see front matter © 2013 Published by Elsevier B.V. http://dx.doi.org/10.1016/j.clay.2013.12.022

Zhuang et al., 2010). Studies show that the contact angle between water and nanocomposite membranes tends to decrease with increasing clay content (Anadao et al., 2010). A systematic study on the nature of clay–water capillary meniscus has not been done in the past. The behavior of clays is controlled at the particle level by the van der Waals (Anandarajah and Chen, 1997), double-layer due to electrostatic interactions (Anandarajah and Chen, 1994; Anandarajah and Lu, 1992), and capillary interactions (Amarasinghe and Anandarajah, 2011; Anandarajah and Amarasinghe, 2012). Particle-level studies (Anandarajah, 1994) require quantification of the interparticle forces that arise from these interactions. While classical equations such as the Young–Laplace equation (Laplace, 1806; Young, 1805) may be used to calculate the capillary forces on clay particles, (a) they require knowledge of contact angles and meniscus curvatures and (b) their validity has not been verified for clay–water systems. The present study uses molecular dynamics (MD) to develop the required understanding, and to verify the validity of the Young–Laplace equation. The capillary phenomenon that occurs in the interparticle space (not interlayer space) is of interest.

1.1. The Young–Laplace equation Capillary meniscus and wetting are primarily characterized by the contact angle θ. Contact angle is the angle at which the liquid–air interface meets the solid surface. The contact angle is a measure of the wettability of the fluid on the solid surface, with θ = 00 representing perfect wettability and θ = 1800 representing perfect non-wettability. For example, θ ≈ 00 for glass–water–air interface and θ ≈ 1400 for glass–mercury–air interface.

P.M. Amarasinghe et al. / Applied Clay Science 88–89 (2014) 170–177

Consider the two-dimensional air–water meniscus between two parallel clay particles as shown in Fig. 1. With the assumption that the meniscus may be approximated by a circular arc of radius r, the Young–Laplace equation for the pressure difference between air and water is pair −pwater ¼

T r

ð1Þ

where pwater is the water pressure inside the meniscus, pair is the air pressure outside the meniscus, T is the air–water surface tension and r is the meniscus radius of the curvature. A positive radius of curvature (Fig. 1) indicates that pwater b pair. The water pressure is negative when the air pressure is zero. The top clay particle is subjected to the following surface tension force in the z − direction: F Tz ¼ −TLx sinθ

ð2Þ

where Lx is the length of the clay layer in the x − direction. Assuming that the pressure is uniform in water, it follows that the clay particles are also subjected to a force Fpz due to the water pressure acting on its side facing the water. Referring to Fig. 1, Fpz is given by

171

There are several public domain and commercial computer codes for performing MD, each with specific pre- and post-processing capabilities. The present study is performed using the computer code named NAMD (Kalé et al., 1999), along with the interactive visual molecular dynamics software named VMD (Humphrey et al., 1996). The theoretical details of NAMD may be found in Phillips et al. (2005). One of the force fields supported by NAMD is the CHARMm (Brooks et al., 1983) field. In it, the total potential energy Utot is expressed as a sum of the bonded energy Ub and the non-bonded energy Unb (i.e., Utot = Ub + Unb). The present study employs the following expression for Unb: U nb ¼

X

X

vanderWaals i≠ j

" 4εij

σ ij r ij

!12 −

σ ij r ij

!6 # þ

X

X qi q j

electrostatic i≠ j

r ij

ð6Þ

where rij is the distance between atoms i and j, εij and σij are constants governing the van der Waals energy, and qi and qj are the electrostatic charges on atoms i and j respectively. Hence NAMD employs the standard Lennard–Jones 6–12 potential for the van der Waals part. The parameters between two dissimilar atoms i and j are computed as follows

ð3Þ

εij ¼

pffiffiffiffiffiffiffiffi εi ε j

ð7aÞ

where Ly is the length of the clay particle that is in contact with water. Combining Eqs. (1–3), the total force Fz acting on the top clay particle due to capillary water then is

σ ij ¼

i 1h σi þ σ j 2

ð7bÞ

F pz ¼ pwater Lx Ly


F z ¼ F Tz þ F pz ¼ −TLx sinθ þ pwater Lx Ly ¼ −TLx

Ly sinθ þ r

 ð4Þ

Note that the sign of the second term will be opposite if the radius of curvature is negative (as in glass–mercury–air meniscus). When the meniscus can indeed be approximated by a circular arc as done here, the radius is related to (half) the particle separation as r = d/cos θ. Eq. 4 then may be written as
 Ly F z ¼ −TLx sinθ þ cosθ d

ð5Þ

It follows from Eq. 5 that the quantification of capillary forces on clay particles in a partially saturated clay–water system requires knowledge of contact angles and meniscus curvatures.

where (εi,σi) are parameters associated with atom i and (εj,σj) are parameters associated with atom j. The molecules involved in the present study are treated as rigid, and hence the bonded energy equations are irrelevant. 2.2. Atomic structures of clay and water The clay used in this study is pyrophyllite, which has the basic structure of smectite type clay mineral. Pyrophyllite has the chemical formula of Al2Si4O10(OH)2 (Van Olphen, 1977). Pyrophyllite is a 2:1 clay mineral. The dimension of a unit cell (which consists of 40 atoms) in the x − y − z Cartesian coordinate system is 5.28 Å × 9.14 Å × 6.56 Å. A typical clay layer consists of several unit cells that are repeated both in the x − and y − directions (Fig. 1). In nature several of these layers are stacked up in the z − direction at some basal spacing. To consider the effect of particle thickness, a single-layer particle and a two-layer particle were analyzed in the present study.

2. Details of molecular dynamics 2.3. Initial configurations and potential energy parameters 2.1. Potential energy functions and simulation software The fundamentals of molecular dynamics are well established; the details are available in many books (e.g., Allen and Tildesley, 1987).

Fig. 1. Schematic of a meniscus between two parallel clay particles.

The initial coordinates of the atoms in the pyrophyllite unit cell were obtained from Skipper et al. (1995). The values of the van der Waals parameters (εij and σij) and the electrostatic charges (qi) used in the study originated from the work of Teppen et al. (1997). These parameters were later cast within the framework of the CHARMm force field, where the functional form of the van der Waals equation is different from those of Teppen et al. (1997) by Katti et al. (2005). These parameters are listed in Table 1. Teppen et al. (1997) have verified the validity of the parameter values by comparing the experimentally measured d(001) spacings with those calculated by molecular dynamics. Thompson et al. (1993) studied the solid/fluid capillary interaction with the flexibility of the solid atoms modeled using springs. The study indicated that stress build up in the fluid was never large enough to produce substantial deformations of the walls. The focus of the current paper is on the calculation of capillary forces between two parallel particles at fixed separations. This requires that the particles be held at fixed locations. In the current study, it was assumed that the change in the initial coordinates of the atoms of the clay layers during the simulation arising from the clay–water interaction is small. On this basis,

172

P.M. Amarasinghe et al. / Applied Clay Science 88–89 (2014) 170–177

Table 1 Van der Waals parameters for clay and water. Atom

εi (kcal/mol·Å)

σi (Å)

qi (e)

Al Si O(interior-1) O(interior-2) O(surface) H(clay) H(water) O(water)

0.150 0.001 6.0 6.0 1.0 0.0001 0.046 0.152

6.30 7.40 2.80 2.80 3.0 2.40 0.44 3.53

1.68 1.40 −0.96 −0.91 −0.70 0.40 0.417 −0.834

the atoms of the clay layers were kept fixed during the simulations. The exact final locations of the fluid atoms near the surface of the particles are likely to depend on whether the clay particle flexibility is modeled or not. In the current study, the meniscus curvature and contact angle were established by focusing on the overall curvature, as was done in previous studies on capillary phenomena (e.g., Martic et al., 2002). The water molecule is modeled using the TIP3P model (Jorgensen et al., 1983), where the molecule is treated as rigid. The most important interaction relevant to the present study is the non-bonded interaction. The electrostatic charges within a unit cell of clay layer sum up to zero. The specific TIP3P model used in the present study is the one employed by CHARMm (Brooks et al. 1983), where the van der Walls parameters are specified separately for the oxygen and hydrogen atoms. The electrostatic charges within each water molecule also sum up to zero. 3. Simulation details The schematic illustration of the construction of the model used in the MD simulations is presented in Fig. 2. Two clay particles are placed parallel to each other keeping a selected spacing (lz) between them. Since the intent is to study the capillary effects in the interparticle space, the value of lz is chosen to be greater than the well-hydrated interlayer spacing of 15.23 Å (Amarasinghe et al., 2008) to avoid any possible interlayer interactions. A rectangular water box of a certain length is placed between the particles. Since, the number of nitrogen and oxygen molecules, which are the most abundant components in the air, occupied in the free space of the simulated model is extremely small (approximately 0.0004 nitrogen molecules and 0.000084 oxygen molecules), air molecules are not included in the simulation. It was assumed that the subsequent change in the interatomic forces is negligible. Periodic boundary conditions are applied in all three directions. Continuous nature of the clay particles and the water body in x − direction is mimicked by keeping the width of the clay particle (lx) exactly the same as the cell size (lx′) in x − direction. The periodic cell size is at least four (van der Waals) cut-off distances larger than the model size in the

• Set 1: Study of the effect of particle thickness • Set 2: Study of the effect of particle spacing • Set 3: Study of the effect of the strength of the clay–water van der Waals attractions

l y + Δ ly ly

clay layer lz + Δlz

lx

water body

lz

z x

y

y − and z − directions (i.e., Δly ≥ 4 × cut off distance and Δlz ≥ 4 × cut off distance) to avoid any influence from the atoms in the neighboring cells. The following procedure was followed in the simulations. First, an independent water box of size 21.12 Å × 270 Å × 270 Å was equilibrated in the isobaric–isothermal ensemble (constant number, pressure, and temperature or NPT) with zero pressure and a temperature of 298 K. Next, rectangular shaped water bodies were cut from the equilibrated water box and placed in between the two clay particles (Fig. 2). The thickness of the water body in the z − direction is chosen so as to fit the water body in the space between the clay particles. To minimize the instability at the beginning of the simulation, a distance of 1–1.5 Å was kept in between the water body and the clay layer. In most cases reported here, the length of the water body was chosen to be 160 Å, except in one case (the large model to be described later), where the length of the water body was chosen to be 200 Å. The number of atoms in the body of water varies with the model used for the study. As described later, three models were used (small, medium and large). The numbers of atoms in the body of water in these models were 12,733, 26,094 and 66,048 respectively. A time step of 1 fs (10−15 s) was used in all MD simulations. Preliminary studies indicated that the time step of 1 fs was optimal for the systems analyzed. A cut-off distance larger than 8 Å is generally considered as a sufficient distance for truncating the van der Waals forces (Cygan, et al., 2004; van Duin and Larter, 2001; Zaidan et al., 2003). Therefore, a cut-off distance of 10 Å was used in this study. The procedure known as the “Ewald sum” was used to calculate the long-range electrostatic forces accurately. The MD simulations were then carried out in the canonical ensemble (constant number, volume, and temperature or NVT) with a temperature of 298 K. With respect to the x − y − z − coordinate system, the size of a clay layer used in the analysis was 21.12 Å × 411.3 Å × 6.56 Å. The layer involves four unit cells in the x − direction (which satisfies the minimum image convention that the cell dimension be greater than twice the cut-off distance), and 45 unit cells in the y − direction. The total number of atoms in a layer was 14,400. In the study employing two layers to represent a particle, a basal spacing of 10 Å (which corresponds to a surface-to-surface spacing of 3.44 Å) was used. All simulations were run until the system attained equilibrium. When the changes in the average values of system energy and capillary force became negligible, the system was considered to have attained equilibrium. To establish equilibrium based on the chosen convergence criteria, the simulation times varied between 2 and 3 ns, depending on the system parameters such as the particle spacing and van der Waals parameters. Three sets of simulations were performed as follows:

lx' = lx

Fig. 2. Placement of the clay–water system in the simulation cell for molecular dynamics study.

The focus of the study reported in this paper was on capillary phenomenon. The intent of Set 1 was to determine the minimum number of clay layers to be used to model a particle so as to study the capillary phenomenon. The study also yields results on effects of particle thickness on capillary phenomenon. Although, the particles are in fixed positions, the non-bonded interatomic forces could affect the shape and the position of the meniscus with increased particle thickness. However, considering the cut-off value of 10 Å used in the study, it is reasonable to assume that the contribution from the van der Walls forces due to increased particle thickness is negligible. Nevertheless, since, the electrostatic interactions are long range, it is reasonable to expect changes in the shape and the position of the meniscus with increased particle thickness. Two simulations were performed, one employing a single layer to model each clay particle, and another employing two layers to

P.M. Amarasinghe et al. / Applied Clay Science 88–89 (2014) 170–177

Slice 10 Slice 9 Slice 8 Slice 7 Slice 6 Slice 5 Slice 4 Slice 3 Slice 2 Slice 1

z y

Slab 4 Slab 3 x Slab 2 Slab 1

Fig. 3. Schematic of slicing the water box for determination of the capillary meniscus.

model each clay particle. A particle spacing of 40 Å was used in this series. In clays found in nature, the inter-particle spacing varies widely. The intent of Set 2 was to study the effect of particle spacing on capillary phenomenon. Three simulations were performed. The values of the spacing used were 40 Å (small model), 80 Å (medium model) and 160 Å (large model). The total number of atoms involved in the systems were as follows: small model: 27,183; medium model: 40,494 and large model: 80,448. The intent of Set 3 was to study the consistency of MD results. The wetting/non-wetting phenomena depend on the strength of the solid/ fluid attraction. To study this, two additional simulations were performed using the model that uses a single layer to model the clay particle and an inter-particle spacing of 40 Å. The LJ potential parameter εi (Table 1) was varied for the oxygen atoms in the clay particles. In one run, εi was increased (from those listed in Table 1, which represents real pyrophyllite) by a factor of 10 and in another, it was decreased by a factor of 10. 3.1. Analysis of results

173

meniscus. Various methods have been used in the past. De Ruijter et al. (1999) used the variation of density in the radial direction as a guide to determine the meniscus. The basis is that the molecular density drops from the equilibrium fluid density to zero from a point in the fluid to air across the meniscus. For a two dimensional configuration such as the ones involved in the current study, one could graphically draw a fitting curve directly on the final configuration. However, although the degree of fluctuation is small compared to quantities such as the system energy, the atomic positions do fluctuate with time. Hence the positions must be time averaged (as is done for the system energy). A preliminary study indicates that the differences between the meniscus determined from the above two procedures are negligible. The last of the two methods is chosen for the present study. Only the last few (say, n) frames were used in averaging. Let us consider frame number one (i.e., first of the last n frames). The objective is to plot a two-dimensional curve of meniscus in the y − z plane. However, the configurations were three dimensional. Hence a spatial averaging needs to be done in the x − direction. With this in mind, the body of water was divided into nx slabs of equal thickness (5.28 Å) in the x − direction. Each slab was then divided into nz equal slices in the z − direction. An example of the divisions with nx = 4 and nz = 10 is shown in Fig. 3. Thus, at a given value of z, there are nx slabs in the x − direction. At this z, the y − coordinate of the outer most atom was found for each of the nx slabs, and averaged to find a single value for the y − coordinate. This was repeated for each of the nz slices. Hence the two-dimensional meniscus in the y − z plane was defined for frame number one. This procedure was repeated for the other n − 1 frames. There was n number of y − coordinates for each value of z. These n numbers of y − coordinates were averaged to find the time average of the y − coordinate for the chosen slice in the z − direction. This was repeated for all nz slices in the z − direction, obtaining the time-averaged, two-dimensional menisci. The y − coordinates determined by the above procedure were then plotted with the z − coordinate of the midpoint of each of the nz slices in the z − direction. The average points representing the meniscus were then fitted with a circular arc, and the contact angle was then obtained by drawing a tangent to the arc at the solid/fluid contact. The procedure was repeated for different values of n in the range of 5 to 100. The results were found to be independent of n provided n≥ 10.

The MD results have been processed to obtain the following: a. The configurations at various times b. The variations of system energy and force vector on a clay particle arising from the clay–water interaction c. The time averages of the system energy and force vector referred to in (b) above d. The time averages of the contact angles and meniscus radius of curvatures The results have been output every 500 time steps (i.e., every 0.0005 ns). The spatial coordinates of atoms, system energy and components of force vectors, among other things, were then available from the start for every 0.0005 ns. The forces of interest are the non-bonded forces between the group of atoms in a given clay particle (top or bottom) and the group of atoms in the body of water. Since all variables fluctuate with time in an MD analysis, an effective method of determining the equilibrium state was required. The time averages of the energy and forces were used to determine this equilibrium point. At equilibrium, the system reaches a “steady state” in terms of these averages; i.e., the averages do not vary with time. Since our primary interest was forces, it was necessary to continue the simulations until not only the energy averages, but also the force averages reach constant values. In order to calculate the capillary forces using the Young–Laplace method, the contact angles and the meniscus radius of curvatures of the air–water interface were needed. Since the contact angle depends on the shape of the meniscus, one has to first accurately locate the

4. Results and discussion 4.1. Determination of equilibrium state Fig. 4 presents the initial and the final configurations for the models with a particle spacing of 40 Å (small model) and regular LJ parameters (Table 1). The model in Fig. 4 uses two layers for the clay particles. The final configurations correspond to a simulation time of 2 ns. Let us investigate the state of equilibrium.

40 Å Initial

final Fig. 4. Initial and final configurations of the two-layer model (2d = 40 Å) indicating the capillary meniscus.

174

P.M. Amarasinghe et al. / Applied Clay Science 88–89 (2014) 170–177

a) Energy (t = 0 – 1.5 ns)

b) Time-averaged energy

Time (ns) 0.5

1

Time (ns) 1.5

-1.8605

Energy x 106 kcal/mol

Energy x 106 kcal/mol

0

-1.8615

-1.8625

-1.861

0.5

1

1.5

-1.862

d) Time-averaged force Fz -Average (kcal/mol.Å)

c) Force (t = 0-3 ns)

Fz (kcal/mol.Å)

0

50

-50

20 -5 -30 -55 -80

-150 0

1

2

Time (ns)

3

0

1

2

3

Time (ns)

Fig. 5. Time variations of energy and force. (a) Actual system energy, (b) average system energy, (c) actual system force, and (d) average system force.

The variation of system energy E with the simulation time t is shown in Fig. 5a. The system energy E jumped as soon as MD simulations began, and gradually decayed thereafter. The variation of E repeated on a regular pattern (Fig. 5a), suggesting that the system has reached equilibrium. However, the variation of the time average of E with time (Fig. 5b) was the one that confirmed that the system has reached equilibrium. The variation of time-averaged E has ceased, indicating that the system has reached a steady state in terms of the time-averaged E. Hence, based on the variation of E, it appeared that a simulation time of about 1 ns was adequate for the system to reach equilibrium. The variation of force Fz and on the top clay particle with time is shown in Fig. 5c. The variation of its time average with time is shown in Fig. 5d. The magnitude of the force fluctuated significantly with time. However, the time-averaged Fz reached a fairly small steadystate value. The steady-state values of Fx, Fy and Fz on the top and the bottom clay particles are listed in Table 2. It may be observed from Table 2 that (1) the steady-state values Fx and Fy are almost zero (relative to the magnitudes of Fz) and (2) the absolute values of the magnitudes of Fz on the top and bottom particles are almost the same. Noting that it takes 1 ns for the time-averaged system energy E to reach a steady state, whereas it takes about 3 ns for the time-averaged

forces to reach steady state, it is clear that the specific quantity of interest (forces, in the present case) and not necessarily the system energy E that must be tracked for terminating the simulation. Due to symmetry, the steady-state values of Fx and Fy were expected to be zero. The MD results confirmed this. The model (Fig. 4) was symmetric about the centerline between the particles and hence the magnitude of Fz on the top particle was expected to be negative of the magnitude of Fz on the bottom particle. Again the MD results were in agreement with this.

Table 2 Comparison of capillary forces on clay particles calculated using molecular dynamics (MD) and Young–Laplace equation. Force (× 10−10) N

By MD simulation Top particle

Bottom particle

Average

Fx

0.020

0.1

Fy

0.08

0.04

Fz

8.21

8.51

(0.02 + 0.1) / 2 = 0.06 (0.08 + 0.04) / 2 = 0.06 (8.21 + 8.51) / 2 = 8.37

By Young–Laplace equation 0 0 9.87

Fig. 6. Clay–water meniscus from MD analysis (× — averaged from last 10 frames and ■ — averaged from last 90 frames) and its approximation by a circular arc (solid line). Meniscus radius is 34.8 Å and contact angle is 55°.

P.M. Amarasinghe et al. / Applied Clay Science 88–89 (2014) 170–177

175

40 Å

40 Å a) Small

Two-layer model 40 Å z

y

80 Å Single-layer model b) Medium

Fig. 7. Comparison of the final configurations of the two-layer and singe-layer models.

4.2. Air–water meniscus Using the procedure described earlier, the air–water meniscus was determined for the two-layer particle model (Fig. 4). The meniscus at the left end was considered. The results (obtained with nx = 4 and nz = 20) are presented in Fig. 6. Also shown in the figure is a circular arc used to fit the meniscus data. The meniscus determined using the last 10 and 90 frames are shown in the figure. The same circular arc is seen to fit the MD results using both the last 10 and 90 frames. This indicates that the last 10 frames are adequate to determine the average meniscus from the results of MD analyses. The radius of the curvature of the meniscus was found to be 34.8 Å. The contact angle was about 55°. This value was very close to the experimentally measured contact angle (between 50° and 60°) for talc (Douillard, et al., 2002), which is another 2:1 phyllosilicate similar to pyrophyllite. The measurements were made by optical observations for water in contact with a talc plate. The receding and advancing contact angles were measured to be 50° and 60° respectively. 4.3. Verification of the validity of the Young–Laplace equation The magnitude of Fz on one of the clay particles was calculated using Eq. 5 for the two-layer particle model shown in Fig. 5. The following parameters were used: T = 0.075 N/m (well-known value for air–water interface), θ = 55ο, Lx = 21.12 Å and Ly = 160 Å. The value calculated by the Young–Laplace equation (9.87 × 10−10 N) compared very well with the average value (time averaged steady-state value of Fz) for the top and bottom particles (8.37 × 10− 10 N) from the MD analysis. The results are listed in Table 2. The results verified the validity of the Young–Laplace equation for the clay–water–air system (as would the additional results presented in the subsequent sections). 4.4. The effect of particle thickness on capillary phenomena The comparison of the final configurations for the two-layer particle model and the single-layer particle model (Fig. 7) suggested that the menisci for the two models are approximately the same. The results presented in Table 3 confirmed this in a quantitative manner. Presented in Table 3 are the time-averaged steady-state values of the components of the forces (average for the top and bottom particles). The results Table 3 Capillary forces for single-layer and two-layer particle models using molecular dynamics (MD) and Young–Laplace equation. Force (× 10−10) N Fx Fy Fz

By MD simulation

By Young–Laplace equation

Single layer

Two layer

Single layer

Two layer

0.07 0.06 9.65

0.06 0.06 8.37

0 0 9.87

0 0 9.87

160 Å

z

c) Large

y

Fig. 8. Effect of particle spacing (40 Å, 80 Å and 160 Å) on air–water meniscus.

indicated that the use of a single layer to model the clay particle was adequate and that the particle thickness has very little influence of the capillary phenomena. On this basis, to reduce the computation time, only the single-layer models were used in the simulations discussed in the following sections.

4.5. The effect of particle spacing on capillary phenomena Simulations were performed on three models (small, medium and large), each with a different particle spacing, to investigate the effects of the particle spacing on capillary phenomena. The final configurations (after a simulation time of 2 ns) are shown in Fig. 8. It may be observed that the contact angles remained the same (550) for the three models, and the radius of meniscus curvature increased as the particle spacing increased. Using the contact angle in the Young–Laplace equation (Eq. 5), the magnitudes of the components of the force vector on the particles were calculated and compared with those from MD in Table 4. Note that d denotes half the spacing between the particles (Fig. 1). The agreements between the MD results and those calculated by the Young–Laplace equation were very good. The force in the z − direction Fz decreased as the spacing increased. This was due to the fact that as the spacing increased, the radius of curvature increased, which

Table 4 Capillary forces on pyrophyllite at different particle spacings calculated from molecular dynamic (MD) analysis and Young–Laplace equation. d, θ0, Ly, and Fz are the half the particle spacing, contact angle, length of the water body and z-direction capillary force respectively. 2d (Å)

θ0

Ly — final (Å)

Fz — average (× 10−10)N MD simulation

Young–Laplace

40 80 160

55 55 55

160 160 200

9.65 6.35 5.20

9.87 6.20 5.00

176

P.M. Amarasinghe et al. / Applied Clay Science 88–89 (2014) 170–177

Table 5 Capillary forces on pyrophyllite calculated from molecular dynamic (MD) analysis using different values of van der Waals parameter εi. θ0, Ly, and Fz are the contact angle, length of the water body and z-direction capillary force respectively. Simulation type

Medium LJ High LJ Low LJ

θ0

εi (kJ/mol) O(interior)

O(surface)

6.0 60.0 0.60

1.0 10.0 0.10

55 0 127

Ly — final (Å)

160 160 160

Fz — average (× 10−10) N MD simulation

Young– Laplace

8.36 14.06 −3.61

9.08 12.67 −4.66

decreased the capillary pressure pwater (Eq. 1) and in turn the force associated with the capillary pressure Fpz (Eq. 3). 4.6. The effect of the strength of van der Waals attraction on capillary phenomena The εi values (LJ parameter) of pyrophyllite oxygen atoms are far greater than those of the other atoms (Table 1). Hence, the pyrophyllite oxygens contribute much more to net clay–water attraction than the other atoms and only the oxygen εi values were changed in this study. The εi value for oxygen was increased by a factor of 10 in one simulation and decreased by a factor of 10 in another simulation for a single-layer particle model with a 40 Å spacing (Table 5). The final configurations (after a simulation time of 2 ns) are compared in Fig. 9. The MD simulations produce the expected results: (a) the increased value of εi, which effectively increased the clay–water van der Waals attraction, decreased the contact angle to almost zero, effectively manifesting “wetting” phenomena (as in glass–water capillary), and (b) the decreased value of εi, which effectively decreased the clay–water van der Waals attraction, increased the contact angle to more than 90° (to 127°), effectively manifesting “non-wetting” phenomena (as in glass–mercury capillary). The force Fz calculated by the Young–Laplace equation (Eq. 5) compared very well with those directly calculated by the MD analysis (Table 5). 5. Summary and conclusions While the phenomenon of capillarity has been studied in the past for various interfaces, the clay–water–air system has remained neglected. The study described here focused on this. The air–water capillary

a) Medium ε i

b) High ε i

z

y

c) Low ε i

Fig. 9. Effect of the van der Waals parameter εi on air–water meniscus.

menisci that form between two parallel clay particles, modeled with the pyrophyllite mineral, were studied using molecular dynamics (MD). The potential energy parameters were obtained from previous studies. The effects of particle thickness, particle spacing and the strength of the clay–water van der Waals attraction on capillary phenomena were investigated. From the MD results, the menisci were determined and the contact angles were estimated. Using this, along with the known value of the air–water surface tension, in the classical Young–Laplace equation, the capillary forces on the clay particles (i.e., the force on the clay particles arising from the clay–water interaction) were calculated, and compared with those obtained directly from the MD analyses. The results were found to be not only consistent with our current understanding based on previous experimental and MD studies, but were in agreement with the results calculated by the classical Young–Laplace equation. More importantly, the study added to the much needed fundamental understanding of the clay–water–air capillary phenomenon.

Acknowledgments The financial support for the study was provided by grants from the National Science Foundation (CMMI0758268 and its REU supplements and CMMI1030570). The supports of Drs. Richard J. Fragaszy and John L. Daniels are acknowledged.

References Allen, M.P., Tildesley, D.J., 1987. Computer Simulation of Liquids. Clarendon Press, Oxford (385 pp.). Amarasinghe, P.M., Anandarajah, A., 2011. Influence of fabric variables on clay water–air capillary meniscus. Can. Geotech. J. 48, 1–8. Amarasinghe, P.M., Katti, K.S., Katti, D.R., 2008. Molecular hydraulic properties of montmorillonite: a polarized Fourier transform infrared spectroscopic study. Appl. Spectrosc. 62, 1303–1313. Anadao, P., Sato, L.F., Wiebeck, H., Valenzuela-Diaz, F.R., 2010. Montmorillonite as a component of polysulfone nanocomposite membranes. Appl. Clay Sci. 48, 127–132. Anandarajah, A., 1994. Discrete element method for simulating behavior of cohesive soils. J. Geotech. Eng. Div. ASCE 120, 1593–1615. Anandarajah, A., Chen, J., 1994. Double-layer repulsive force between two inclined platy particles according to Gouy–Chapman theory. J. Colloid Interface Sci. 168, 111–117. Anandarajah, A., Chen, J., 1997. Van der Waals attractive force between clay particles in water and contaminant. Soils Found. Jap. Soc. Soil Mech. Found. Eng. 37, 27–37. Anandarajah, A., Lu, N., 1992. Numerical study of the electrical double-layer repulsion between non-parallel clay particles of finite length. Int. J. Numer. Anal. Methods Geomech. 15, 683–703. Anandarajah, A., Amarasinghe, P.M., 2012. Microstructural investigation of soil suction and hysteresis of fine-grained soils. J. Geotech. Geoenviron. Eng. ASCE 138, 38–46. Bitsanis, I., Magda, J., Tirrell, M., Davis, H.J., 1987. Molecular dynamics of flow in micropores. Chem. Phys. 87, 1733–1750. Brooks, B.R., Bruccoleri, R.E., Olafson, B.D., States, D.J., Swaminathan, S., Karplus, M.J., 1983. CHARMm: a program for macromolecular energy, minimization, and dynamics calculations. J. Comput. Chem. 4, 187–217. Cipriano, B.H., Raghavan, S.R., McGuiggan, P.M., 2005. Surface tension and contact angle measurements of a hexadecyl imidazolium surfactant adsorbed on a clay surface. Colloids Surf. A Physicochem. Eng. Asp. 262, 8–13. Cygan, R.T., Liang, J., Kalinichev, A.G., 2004. Molecular models of hydroxide, oxyhydroxide, and clay phases and the development of a general force field. J. Phys. Chem. B 108, 1255–1266. De Ruijter, M.J., Blake, T.D., De Coninck, J., 1999. Dynamic wetting studied by molecular modeling simulations of droplet spreading. Langmuir 15, 7836–7847. Docoslis, D., Giese, R.F., van Oss, C.J., 2000. Influence of the water–air interface on the apparent surface tension of aqueous solutions of hydrophilic solutes. Colloids Surf. B: Biointerfaces 19, 147–162. Douillard, J.M., Zajac, J., Malandrini, H., Clauss, F., 2002. Contact angle and film pressure: study of a talc surface. J. Colloid Interface Sci. 255, 341–351. Heinbuch, U., Fischer, J., 1989. Liquid flow in pores: slip, no-slip, or multilayer sticking. Phys. Rev. A 40, 1144–1146. Humphrey, W., Dalke, A., Schulten, K., 1996. VMD — visual molecular dynamics. J. Mol. Graph. 14, 33–38. Ichikawa, Y., Kawamura, K., Fujii, N., Kitayama, K., 2004. Microstructure and micro/macrodiffusion behavior of tritium in bentonite. Appl. Clay Sci. 26, 75–90. Jorgensen, W.L., Chandrasekhar, J., Madura, J.D., 1983. Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 79, 926–935. Kalé, L., Skeel, R., Bhandarkar, M., Brunner, R., Gursoy, A., Krawetz, N., Phillips, J., Shinozaki, A., Varadarajan, K., Schulten, K., 1999. NAMD2: greater scalability for parallel molecular dynamics. J. Comput. Phys. 151, 283–312.

P.M. Amarasinghe et al. / Applied Clay Science 88–89 (2014) 170–177 Katti, D.R., Schmidt, S.R., Ghosh, P., Katti, K.S., 2005. Modeling response of pyrophyllite clay interlayer to applied stress using steered molecular dynamics. Clay Clay Miner. 53, 171–178. Koplik, J., Banavar, J., Willemsen, J., 1989. Molecular dynamics of fluid flow at solid surfaces. Phys. Fluids A 1, 781–794. Laplace, P.S., 1806. Supplement to book 10 of Mecanique Celeste. Lee, J.Y., Lee, S.H., Kim, S.W., 2000. Surface tension of silane treated natural zeolite. Materi. Chem. Phys. 63, 251–255. Martic, G., Gentner, F., Seveno, D., Coulon, D., De Coninck, J., 2002. A molecular dynamics simulation of capillary imbibitions. Langmuir 18, 7971–7976. Nakano, M., Kawamura, K., Ichikawa, Y., 2003. Local structural information of Cs in smectite hydrates by means of an EXAFS study and molecular dynamics simulations. Appl. Clay Sci. 23, 15–23. Phillips, J.C., Braun, R., Wang, W., Gumbart, J., Tajkhorshid, E., Villa, E., Chipot, C., Skeel, R.D., Kalé, L., Schulten, L., 2005. Scalable molecular dynamics with NAMD. J. Comput. Chem. 26, 1781–1802. Skipper, N.T., Chang, F.R., Sposito, G., 1995. Monte Carlo simulations of interlayer molecular structure in swelling clay minerals. 1. Methodology. Clay Clay Miner. 43, 285–293. Sokołowski, S., 1991. Capillary condensation and molecular dynamics simulations of flow in narrow pores. Phys. Rev. A 44, 3732–3741.

177

Teppen, B.J., Rasmussen, K., Bertsch, P.M., Miller, D.M., Schafer, L., 1997. Molecular dynamics modeling of clay minerals. 1. Gibbsite, kaolinite, pyrophylite and beidellite. J. Phys. Chem. B 101, 1579–1587. Thompson, P.A., Brinckerhoff, W.B., Robbins, M.O., 1993. Microscopic studies of static and dynamic contact angles. J. Adhesion Sci. Technol. 7, 535–554. Thompson, P.A., Robbins, M.O., 1989. Simulations of contact-line motion: slip and the dynamic contact angle. Phys. Rev. Lett. 63, 766–769. Trozzi, C., Ciccotti, G., 1984. Stationary nonequilibrium states by molecular dynamics. II. Newton's law. Phys. Rev. A 29, 916–925. Van Duin, A.C.T., Larter, S.R., 2001. Molecular dynamics investigation into the adsorption of organic compounds on kaolinite surfaces. Org. Geochem. 32, 143–150. Van Olphen, H., 1977. An Introduction to Clay Colloid Chemistry. Wiley Interscience, New York. Young, T., 1805. An essay on the cohesion of fluids. Philos. Trans. R. Soc. Lond. 95, 65–87. Zaidan, O.F., Greathouse, J.A., Pabalan, R.T., 2003. Monte Carlo and molecular dynamics simulation of uranyl adsorption on montmorillonite clay. Clay Clay Miner. 4, 372–381. Zhuang, J., Goeppert, N., Tu, C., McCarthy, J., Perfect, E., McKay, L., 2010. Colloid transport with wetting fronts: interactive effects of solution surface tension and ionic strength. Water res. 44, 1270–1278.