Fluid Phase Equilibria 501 (2019) 112253
Contents lists available at ScienceDirect
Fluid Phase Equilibria j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / fl u i d
Effect of alkyl chain length on the wetting behavior of imidazolium based ionic liquids: A molecular dynamics study Sanchari Bhattacharjee, Sandip Khan* Department of Chemical & Biochemical Engineering, Indian Institute of Technology Patna, Patna, 801103, India
a r t i c l e i n f o
a b s t r a c t
Article history: Received 19 January 2019 Received in revised form 19 July 2019 Accepted 23 July 2019 Available online 26 July 2019
Molecular dynamics (MD) simulations are performed to investigate the wetting behavior of imidazolium-based ionic liquids (ILs), which mainly consist of an imidazolium ring with a flexible alkyl chain attachment. The length of the alkyl chain can be easily manipulated to modulate the wetting properties of the ILs. To understand the role of the alkyl chain, we investigate the wetting behavior of two room-temperature ionic liquids (RTILs), 1,3-dimethylimidazolium tetra-fluoroborate [DMIM][BF4] and 1propyl-3-methylimidazolium tetra-fluoroborate [PrMIM][BF4], on a graphite surface. Furthermore, the effects of the droplet size and temperature on the contact angle of the IL droplets are investigated. The wetting behaviors of both IL droplets are characterized using a density profile and orientation order parameter profile along the axis normal to the surface. The length of the alkyl chain is found to be crucial for modulating the wetting properties of the ILs. With an increase in the length of the alkyl chain, the contact angle of the droplet decreases, which is also in line with experimental results [1,2]. The wetting of the IL droplets strongly depends on the temperature and the relative fluid-fluid and fluid-substrate interactions. With the increase in temperature, the contact angle of the IL droplet decreases, resulting in the formation of a single layer film at a high temperature. In addition, the behavior of the contact angle of the IL droplet is analyzed through the surface tension as well as the graphite-IL interfacial tension. Furthermore, considering the anisotropic properties of the imidazolium-based ILs, we calculate the contact angles for cylindrical filaments (2D) on the graphite surface and compare the obtained values with the contact angle obtained from the spherical droplet (3D). The contact angle of the cylindrical droplet is found to have a negligible system size effect as compared to the case with the spherical droplet. © 2019 Elsevier B.V. All rights reserved.
Keywords: Wetting Ionic liquid Graphite Alkyl chain length Spherical droplet Cylindrical droplet
1. Introduction In recent years, ionic liquids (ILs) have attracted considerable attention due to their wide range of applications [3e6]. ILs are salts that comprise solely of cations and anions. The properties of ILs, therefore, can be modulated to suit a specific application by selecting appropriate cations and anions [7,8]. Due to their excellent performance in various applications including in nanoelectromechanical systems, transportation in capillary porous media, micro-lenses, nano-lubrication, and digital micro/nanofluidic devices [9e11], they have attracted significant attention, and the study of their wettability is receiving growing interest [1,2,12e27]. The fluid behavior at the micro/nanoscale near the surface is quite different from that of the bulk, which is the inherent
* Corresponding author. E-mail address:
[email protected] (S. Khan). https://doi.org/10.1016/j.fluid.2019.112253 0378-3812/© 2019 Elsevier B.V. All rights reserved.
part of any micro/nanoscale device [28,29]. Extensive studies on the wetting behavior of water droplets on different surfaces have been carried out [30e34]; however, few works have been performed on the wetting behavior of ILs using molecular simulation due to some inherent difficulties, including slow dynamics, strong coulombic interactions between cations and anions, bulky nature, and heterogeneity in structure [13,15,19e21,24,25,27,35]. Furthermore, adjacent to the solid surface, an IL droplet forms a very stable layer, which retards the spreading of the three-phase contact line, and occasionally leads to meta-stable wetting states. For example, the final contact angle of the droplet may be different if the two initial configurations are extremely different (one is in the complete wetting state and other is in the non-wetting state), as observed by Burt and co-workers [13]. Moreover, the spreading of the threephase contact line is different in the x and y directions due to their structural heterogeneity, as well as the asymmetric crystal structure of the surface, which results in the non-spherical shape of the droplet [19]. Thus, extracting the liquid-vapor interfacial points
2
S. Bhattacharjee, S. Khan / Fluid Phase Equilibria 501 (2019) 112253
is considerably difficult, and erroneous results for the measurement of the contact angle are usually provided. Therefore, some authors [19,21] have reported different contact angles for the x- and y-directions. The wettability of IL also depends on the system size at the nanoscale due to the strong layering of the ILs near the solid surfaces [13,19,20], and the length of the three-phase contact line is finite. Thus, understanding the wetting behavior of the ILs on the different surfaces at the molecular level is definitely challenging, as evidenced by the few publications in this area. Therefore, efficient techniques are required to characterize the wetting behavior of ILs. In this regards, a cylindrical drop of IL can be an alternate method of considerably overcoming the system size effect, as the three-phase contact line is infinite and can eliminate the different spreading of x and y dimensions as one of these dimensions is infinite. In this method, the molecular arrangement can be easily examined at the solid-fluid and fluid-fluid interfaces, and it is a crucial parameter for understanding the wetting behavior; moreover, the interfacial points can be traced easily to compute an accurate contact angle for characterizing the wetting behavior. Significant progress has been made on a macroscopic level; however, knowledge of the wetting behavior of ILs on the molecular level is very limited [1,2,12,16e18]. Burt and co-workers [13] examined the wetting behavior of [EMIM][BF4] on graphene sheets at different surface potentials and different temperatures and calculated the contact angles by analyzing the two-dimensional (2D) atomic resolution. The contact angle of the droplet decreases as the temperature and surface interaction potential increases. In addition, they reported an optimum surface interaction parameter of graphite based on the system size, above which the droplet completely spreads on the surface at high temperatures. Most importantly, they observed the metastable wetting state at low temperatures with high surface potentials. Herrara et al. [20] investigated the wetting characteristics using [EMIM][GLY] on graphene, silica, and graphene sheet supported on silica; they found that coating significantly influenced the wetting mechanism. In addition, they reported that the contact angle of the droplet strongly depended on the system size in all cases. Castejon et al. [15] studied both the wetting and tribological properties of various ILs on the surface of silica and found that the contact angles of all the ILs employed were significantly lower than 90 . The nature of the anion also plays an important role in the wetting properties of ILs, as observed by Guan et al. [19]. They found that the contact angles of the IL droplets are in the sequence of [EMIM][Cl] >[EMIM] [PF6] >[EMIM][BF4], which correlates with the viscosities of the ILs. The IL droplet with a relatively low viscosity easily spreads on the surface, thereby resulting in a low contact angle of the IL droplet. Conversely, Malali et al. [21] studied the wetting behavior of [BMIM][PF6] on the surface of titanium dioxide (TiO2) by MD simulation. They observed that the different contact angles in the x and y directions of the IL droplet are due to the anisotropic characteristics of the TiO2 surface. They also observed that the
imidazolium ring of the cations is in parallel with the surface due to the strong surface-IL interaction, which results in the formation of a very stable layer adjacent to the surface that has considerably low mobility. From the above discussion on previous studies, it can be realized that the wetting behaviors of imidazolium-based ILs are very complex mainly due to their structural heterogeneity. The imidazolium-based ILs primarily consist of planar imidazolium rings along with flexible alkyl chains. Therefore, expectedly, the flexibility of the imidazolium cation would play an important role in the wetting behavior. Thus, [DMIM] (primarily consists of planar imidazolium ring) and [PrMIM] (consists of relatively long flexible alkyl chains in addition to the imidazolium ring) are employed here, to understand the influence of the structural integrity on the wetting phenomena. In addition, considering the anisotropic behavior of an IL droplet on the surface along with threedimensional (3D) spherical droplets, we have examined the cylindrical droplet and compared the contact angles for both systems. This paper is organized as follows: Section 2 provides a detailed description of the computational simulation, while Section 3 discusses the simulation results. Subsequently, Section 4 provides concluding remarks. 2. Computational model and methods In this study, two kinds of imidazolium-based ILs, [DMIM][BF4] and [PrMIM][BF4], are investigated to understand the wetting behavior on graphite surface using MD simulation. The molecular structures of the cation and anions are shown in Fig. 1. The interaction between the ILs is expressed by Eq. (1):
2 Xbonds kr;ij 2 Xangels kq;ijk U rij ¼ qij q0;ijk r r0;ij þ ij ijk 2 ij 2 i Xdihedrals X4 Vm;ijkl h 1 þ ð 1Þm cos mfijkl þ ijkl m¼1 2 8 9 2 !12 !6 3 XX< s s qi qj = 1 ij ij 5þ þ 4ε 4 : ij 4pε0 rij ; rij rij i
j>i
(1) whererij is the distance between any pair of atoms i and j; s, the atom size; ε, the interaction energy depth between two atoms; qi and qj , the charges centered on the individual atoms of different IL molecules; kr kq , and Vm , the coefficients for bond stretching, bending, and dihedrals, respectively. q is the angle between three atoms (ijk); 4, the dihedral angle; m, the multiplicity of function, and ε0 ; the dielectric constant. The force field parameters were obtained from the work of Agilio Padua and J.N. Canongia Lopes [36,37]. The LJ parameters for the carbon atoms of graphite sheet were obtained from the work of Jorgensen et al. [38] and Maolin et al. [39], with a size parameter (s) of 0.337 nm and well depth (ε)
Fig. 1. Molecular structure of 1,3-dimethylimidazolium [DMIM], 1-propyl-3-methylimidazolium [PrMIM] and tetra-fluoroborate [BF4].
S. Bhattacharjee, S. Khan / Fluid Phase Equilibria 501 (2019) 112253
of 0.2929 kJ mol1. The interfacial properties of the IL on the graphite surface using this interaction parameter for carbon atoms has been reported by many researchers [40,41]; furthermore, the nanowetting behavior of ILs on graphene sheet using the same interaction parameters was recently reported by Cesar Herrera et al. [20]. We adopted the same interaction parameters to work in line with previous studies. The cross interaction parameters were evaluated using the LorentzBerthelot arithmetic mixing rules [42]. A cutoff distant of 12 Å was used for non-bonded interactions. Initial Configuration: The AB stacking of graphene layers was used to model the graphite surface with an interlayer spacing of 3.35 Å. Three layers of the graphene sheet were applied as the surface to cover the entire potential cutoff range. The graphene sheets were maintained in a frozen state during the simulation. All the three directions of the simulation box were considered as periodic. However, the height of the simulation box was taken large enough to avoid any interaction with the periodic image of the droplet. A number of (150e400) ion pairs of ILs were distributed randomly in a cubic box using the Packmol software, depending on the system of interest [43]. Subsequently, the cubic box of the ILs was placed slightly above the graphite surface (within the distance of 10 Å), as shown in Fig. 2a. The xey dimension of the graphite sheet was varied from (10 10 nm2) to (20 20 nm2) depending on the number of ion pairs. During the equilibration period, the ILs gradually turned into droplets on the surface (Fig. 2b). For the cylindrical droplet, the y-dimension (150 Å) was much longer than the x-dimension (50 Å) (Fig. 2c). Initially, the cubic box of the ILs (50 Å) covered the entire length of the x-dimension so that the periodic image of the box in the x-dimension stabilized the liquid filament. During the equilibration period, it slowly turned into liquid filament, as shown in Fig. 2c and d. Simulation Details: All MD simulations were carried out in a canonical ensemble (NVT) using an open-source package LAMMPS software [44] The particleeparticle particleemesh (PPPM) method [45] with an accuracy value of 1.0 e4 was used to calculate the long-range coulombic interaction. The velocity Verlet algorithm was used to determine the trajectories of each atom with an integration time of 1.0 fs. A NoseeHoover [46] thermostat with damping coefficients of 0.1, was used to set the desired temperature in the respective systems, and the simulation temperature ranged from 300 to 400 K for the studied ILs. Initially, we equilibrated the system using a canonical ensemble (NVT) at T ¼ 300 K for 30e35 ns.
3
Once we determined the constant contact angle of the droplet for the system of interest over some period of time (typically 4e5 ns slightly over 30e35 ns), we conducted the simulation once more at a slightly higher temperature (T ¼ 325 K) than previously used and annealed the system to the desired temperature (T ¼ 300 K) for 3e4 ns to avoid any meta-stable configurations during the simulation. A temperature difference of 25 K was chosen, as this difference was enough to enhance the spreading of the droplet significantly. During annealing of the system (for 3e4 ns), initially, the droplet spread increased slightly (i.e., contact angle slightly decreased for T ¼ 325 K) and returned to the initial shape as the temperature of the system gradually approached the desired temperature (i.e., T ¼ 300 K). We repeated this procedure 3e4 times to ensure consistency in the contact angles obtained. Afterward, another 10 ns run was allowed for calculating the properties of the system. The structures were visualized using a molecular graphic program named Visual Molecular Dynamics (VMD) [47]. Estimation of the Contact Angle: The contact angle of a droplet was evaluated by finding a tangent across the vapor-liquid interface at the three-phase contact line. The vapor-liquid interface was traced by 2D binning across the center of the droplet, as described by Ruijter et al. [48] In this approach, first, the center of mass of the droplet was calculated in the xey plane, and the z-axis was fixed through the center of mass of the droplet. Two types of binning were considered to locate the vapor-liquid interfacial points: rectangular binning in the z-direction with a thickness of 1 Å, and radial binning in the xey plane. The pradius ffiffiffiffiffiffiffiffiffiffiffiffi of each bin in the radial direction was estimated as ri ¼ idA=p (for i ¼ 1 to number of radial bin), where dA is 95 Å2, i.e., each bin was of equal volume. The density profile along the radial direction for each z-bin was fitted with Eq. (2) to locate the interfacial points between the liquid and vapor phase.
rðrÞ ¼
1 1 L 2ðr re Þ r þ rV rL rV tanh 2 2 d
(2)
where rL and rV are the liquid and vapor densities, respectively; r, the radial distance from the z-axis for each layer parallel to the xey plane; re ; the interfacial points of the liquid and vapor phases; and d; the interfacial thickness. Finally, the least square method was used to fit the interfacial points obtained from the above procedure into a circle, after which the contact angle was estimated at the
Fig. 2. Snapshot of 250 ion-pairs of [DMIM][BF4] on graphite surface (a) Initial configuration and (b) Final equilibrate shape for spherical droplet (c) Initial configuration and (d) Final equilibrate shape for cylindrical droplet.
4
S. Bhattacharjee, S. Khan / Fluid Phase Equilibria 501 (2019) 112253
three-phase contact line. The reference for evaluating the contact angle for the droplet at the solid-fluid interface was obtained as the atom position of the first layer of the graphite surface, which was fixed during the simulation. However, we omitted the interfacial points up to 10 Å from the surface while fitting the interfacial points into a circle to avoid the density fluctuation adjacent to the surface. Estimation of the surface tension: Three-hundred ion pairs inside a cubic box of edge length 50 Å in the x and y directions and an elongation in the z-axis (four times the base, 200 Å) were considered to study the liquidevacuum interface. A 10 ns equilibration run was employed at 300 K, and another 5 ns was performed to calculate the average surface tension based on the IrvingeKirkwood method [49,50]; three runs were performed for each IL. The interfacial tension between the graphite and liquids was calculated from the measured liquid surface tensions and the contact angle, q, between the liquid and the solid, using the Young equation, gSL ¼ gSV gLV cosðqÞ. We assumed the value from the literature for the surface energy of graphite, gSV ¼ 54:8 mN/m [51]. The surface tension, gLV , was determined by Eq. (3).
gLV ¼
LZ ðhPT ðZÞi hPN ðZÞiÞ 2
(3)
where Lz is the box length; PN ðZÞ ¼ PZZ ðZÞ normal and PT ðZÞ ¼ ðPXX ðZÞ þPYY ðZÞÞ=2 tangential components of the pressure tensor. The wettability of a surface can also be related to the work of adhesion (wd ), which is defined as the change in the interfacial energy per unit area and can be determined from the YoungeDupre equation, as follows [34].
wad ¼ gLV ð1 þ cos qÞ
(4)
3. Results and discussion
Fig. 3. Graphite-IL interaction energy for 250 ion pair of spherical droplet [DMIM] [BF4]. The interaction energy was almost constant during 30e35 ns (labeled as “Constant”, after 30ns). The system was then periodically annealed for next 8e10 ns (marked as “Annealing”, 36e45 ns). Finally, all relevant properties of IL droplet were estimated from 45 ns onwards (labeled as “Production”). Here inset snapshot of droplet “a”, “b” and “c” corresponding to time 5 ns, 36 ns and 45 ns respectively.
constant. The graphite-IL interaction for [DMIM][BF4] was found to be constant after 30e35 ns, after which the system was annealed for another 8e10 ns to determine the final contact angle of the droplet (~71, see inset image in Fig. 3). Furthermore, to ensure no further decrease in the contact angle of the droplet, a separate simulation was performed for the configuration of the droplet with a relatively low contact angle (~50 ), and the final contact angle of the droplet was found to be the same as that obtained in the earlier approach, in which a cube of IL was placed on the surface as starting configuration. Finally, another 10 ns was used for the production of the system to compute the contact angle, density profile, and orientation profile of the IL droplets on the graphite surface.
3.1. Equilibrium of droplet 3.2. Effect of drop size The dynamics of IL were found to be considerably slow due to their complex structures, bulky nature, and strong coulombic interactions between the cations and anions. Therefore, it required a relatively long time to achieve the equilibrium structure near the solid surface depending on the size and nature of the cations and anions present in the IL. Moreover, due to their highly asymmetric structure, cation and anion molecules may be trapped in a local minimized structure, which occasionally takes time to overcome the energy barrier before attaining the equilibrium state of the droplet. In this work, we examined the wettability of [DMIM][BF4] and [PrMIM][BF4] on the graphite surface by observing the contact angles of the droplets. The contact angle was evaluated from the atomic resolution of the droplet, as described in the methodology section. The [DMIM] cation was mainly comprised of a planar aromatic ring, whereas [PrMIM] consisted of a relatively long flexible alkyl chain in addition to the aromatic ring, as shown in Fig. 1. To avoid any meta-stable structure of the system, we periodically annealed the system (after 30e35 ns) at a slightly high temperature (25 K) until the equilibrium structure of the droplet was determined, as discussed in the earlier section. The equilibrium period of the system was traced by observing the graphite-IL interaction to be constant. The graphite-IL interaction was mainly owing to van der Waals forces and the cations were found to have a larger contribution compare to that of the anions, as shown in Fig. 3. The equilibrium period of the droplet was mainly due to the rearrangement of the cations near the surface; noticeably, a relatively long time was required for the cation-graphite surface to be
The contact angle of a droplet generally varies in a linear manner with the increase in the droplet size, as can be found in several studies [13,20]. To investigate the system size effect for our systems, we analyzed the wetting behavior of the [DMIM][BF4] droplet with the different numbers of ion pairs. Fig. 4 shows the
Fig. 4. Variation of contact angle for spherical droplet with the number of ion pairs of [DMIM][BF4]. Error bars represent the standard error (SE) of the contact angle.
S. Bhattacharjee, S. Khan / Fluid Phase Equilibria 501 (2019) 112253
contact angle dependency on the droplet size. The droplet size effect on the contact angle of the droplet was found to be significant when the droplet consisted of less than 200 ion-pairs. The equilibrium contact angle of the droplet was mainly due to the relative interactions between graphite-IL and IL-IL. When the droplet size was considerably small (less than 200 ion-pairs), the graphite-IL interaction dominated the IL-IL interaction, thereby decreasing the contact angle of the droplet. However, at a relatively large droplet size, the system size effect was considerably less, which is consistent with the previous studies [19]. Fig. 5 represents the change in the contact angle with time for a spherical droplet of [DMIM][BF4] for different droplet sizes. Interestingly, the equilibrium periods for all droplet sizes were approximately the same, around 30e35 ns. However, the computational cost was very high for large-sized droplets. Noticeably, the final contact angle for the 250 ion-pair and 390 ion-pair converged to 72±1 . Therefore, the system with 250 ion-pair was taken further to analyze the wetting behavior of the droplet. The contact angle of the droplet for the 250 ion-pair of [PrMIM][BF4] was determined as 62± 1 ; which is considerably lesser than that of [DMIM][BF4] [1,2]. The error bars for the contact angle of the IL droplet were evaluated by taking the standard error (SE) over 1000 configurations at intervals of 4 ps. 3.3. Density profile The characteristic features of the droplet of both the ILs were examined through the density profile along the direction perpendicular to the surface (i.e., along the z-coordinate in our systems). The density profile of the droplet was evaluated for the cylindrical regime of radius 10 Å, maintaining the axis of the cylinder along the z-axis. The axis of the cylinder originated from the center of the mass of the droplet projected on the surface. The cylindrical regime is shown in Fig. 6 as a dotted area inside the droplet. This method was used to extract the actual mass density of the ion pairs within the droplet along the z-coordinate. The atomic co-ordinate of the ion pairs present in the droplet was updated in the cylindrical regime while calculating the density profile. The bin width along the z-axis was taken as 0.5 Å, which is small enough to predict the density profile with reasonable accuracy. The density profile of the equilibrium droplet was calculated for 5 ns at intervals of 0.01 ps. Fig. 7a and b represent the cations, anions, and total mass density profile of the [DMIM][BF4] and [PrMIM][BF4] droplets, respectively.
5
Fig. 6. The Cylindrical regime of radius 10 Å, represented by dotted area, is used to evaluate the density profile of the droplet along the z-axis.
Fig. 7a shows the strong layering of the cations near the graphite surface compared to that of the anions, which have also been reported in various literatures due to the VDW interactions that play a significant role in attracting the ions towards the surface [19e21]. The considerably sharp first peak of the cations appeared at ~3.5 Å, which indicated the strong orientation ordering of the cations near the surface, which will be discussed in the next section. Interestingly, we observed the almost complete separation of the 1st and 2nd layers of the cation [DMIM] indicating the presence of moieties in between the two layers (Fig. 7a, blue line). Conversely, the 1st peak of the anions was quite broad and almost uniformly distributed across the peak (Fig. 7a, red line). The broad peak of the anion was attributed to the tetrahedral structure of the [BF4] anion. Moreover, we observed a relatively weak layering at the vaporliquid interface (away from the surface) due to the specific orientation of the cations at the liquid-vapor interface, which was also observed in previous work [52]. The density profile of [PrMIM][BF4] (shown in Fig. 7b) was similar to that of [DMIM][BF4], although the layering of the cations near the surface was comparatively less. To understand the role of the alkyl chain, we further analyzed the density profile of the [PrMIM] cation close to the surface and plotted the density profile of imidazolium ring and alkyl chain (consists the last two methyl groups) separately (shown in Fig. 8). The alkyl groups were found to be preferentially adsorbed near the surface and in the vapor-liquid interface. Fig. 9 represents the snapshot of the 1st layer near the surface for [DMIM] and [PrMIM] cations. The thickness of the first layer was determined as the distance from the surface to the first minima located in the density profile. We found that the three-phase contact line was not spherical, and was more irregular for the [PrMIM][BF4] IL. This is due to the higher heterogeneity in the structure of [PrMIM] cations in comparison to that in the structure of the [DMIM] cations. Noticeably, most of the cations were parallel to the surface to maximize the interaction with the surface. The alkyl chain for [PrMIM] in the three-phase contact line was found to be oriented towards the vapor phase. This arrangement of the alkyl chain in the [PrMIM] cation was mainly responsible for the enhancement of the wetting behavior of [PrMIM][BF4]. 3.4. Orientation order parameter
Fig. 5. Evolution of contact angle of the droplet with time for the different number of ion pairs of [DMIM][BF4].
The layering effect in the density profile was mainly due to the orientation of cations near the surface. To understand the structure of the cations near the surface, the profile of the average orientation order parameter was examined along the z-coordinate normal to the surface. The average orientation order parameter was evaluated
6
S. Bhattacharjee, S. Khan / Fluid Phase Equilibria 501 (2019) 112253
Fig. 7. Density profile for spherical droplet of 250 ion pair along the height z (a). [DMIM][BF4] (b) [PrMIM][BF4].
to van der Waals forces. Conversely, the positive value of Sz for the second layer implies that the imidazolium ring was slightly tilted towards the z-coordinate, which is also in line with previous work [19,20]. A similar observation was also found for [PrMIM] cations, as shown in Fig. 10b. 3.5. Interaction energy
Fig. 8. Density profile of Imidazolium ring and alkyl chain of 250 ion pairs of spherical droplet [PrMIM][BF4] along height Z.
as follows:
sz ¼
3 cos2 ðqÞ 1 2
(5)
where q is the angle between the z-coordinate normal to the surface and the vector of interest. Fig. 10a depicts the value of the orientation order parameter of the N1_N2 and CR_CW vectors of the [DMIM] cation as a function of the z-coordinate. The negative values of Sz of both vectors in the 1st adsorbed layer suggest that the imidazolium ring tends to reside parallel to the surface to maximize the interaction with the surface, which is mainly owing
The structure of [DMIM][BF4] near the surface was quite different from that of [PrMIM][BF4], as discussed in the earlier section. This is mainly due to the presence of longer flexible alkyl chains attached with an imidazolium ring, for [PrMIM][BF4]. Therefore, expectedly, the interaction with the graphite surface would be different for both ILs. Fig. 11 shows the interaction energy (normalized with the number of IL pairs) between graphite and [DMIM][BF4] with an increase in the system size. Noteworthily, the contact angle of [DMIM][BF4] was almost constant in the range of the droplet size. It can be observed from Fig. 11 that the contribution of the cations to the graphite-IL interaction energy is more than that of the anions. To characterize the wetting behaviors of [DMIM][BF4] and [PrMIM][BF4], we compared the interaction energy of both the ILs with the graphite surface. The surface interaction energy was found to be higher for [PrMIM][BF4] compared to that of [DMIM][BF4], as shown in Fig. 11. This implies that the [PrMIM] cations have a higher affinity towards the graphite surface, which thus corresponds to the decrease in the equilibrium contact angle of [PrMIM][BF4] compared to the case with [DMIM][BF4]. 3.6. Effect of temperature To understand the influence of the temperature on the wetting behaviors of the IL droplet, 250 IL ion pairs of [DMIM][BF4] and [PrMIM][BF4] drops were investigated in the range of 300e400 K
Fig. 9. Snapshot of the first layer of ILs droplet representing the nature of the shape of the base of the droplet and the orientation of cations of 250 Ion-pair for spherical droplet (a) [DMIM] cations (b) [PrMIM] cations (c) Alkyl chain of [PrMIM].
S. Bhattacharjee, S. Khan / Fluid Phase Equilibria 501 (2019) 112253
7
Fig. 10. Orientation order parameter of N1_N2 vector and CR_CW vector of (a) [DMIM] cations (b) [PrMIM] cations for 250 Ion-pair of spherical droplet along the Z-direction.
Fig. 11. Graphite-IL interaction energy per ion pair for [DMIM][BF4] and [PrMIM][BF4].
(shown in Fig. 12). We adopted the equilibrated drop of the 250 ion pairs of [DMIM][BF4] and [PrMIM][BF4] at 300 K as the initial configuration for the simulation for different temperatures at intervals of 25 K. Fig. 12 shows the contact angle variation of the IL droplet with an increase in temperature in the range of 300e350 K. In general, the contact angle of the droplet decreases with an
increase in temperature for both ILs. For each temperature, the contact angle of the [DMIM][BF4,] droplet was higher than that of [PrMIM][BF4]. Furthermore, as the temperature increased, the IL droplet spread increased and formed a thin film at ~400 K. A similar trend of the contact angle with increasing temperature was also observed for other IL droplets [19]. This may be due to the decrease in the surface tension with increasing temperature, as reported by Xu et al. [53] They reported the value of surface tension for various imidazolium-based ILs, [CnMIM][BF4] (n ¼ 2,3,4,5,6), in the temperature range of (298.15e338.15 K) and observed that the surface tension of the ILs decreased as the temperature and chain length increased. A similar observation was found for other ILs in literature [54,55]. Therefore, the decrease in the contact angle of the studied IL droplets on graphite surface was observed as the temperature and chain length increased (Fig. 12). However, the solid-fluid interfacial tension also significantly influences the contact angle of the droplet, which is discussed in the next section. The density profile along the z-coordinate through the center of the droplet was evaluated to examine the effect of temperature (shown in Fig. 13) on the structure of the IL droplet. The number of layers normal to the surface decreased with increasing temperature, indicating the spreading of the droplet on the surface, which implies a decrease in the contact angle of the droplet, thereby forming a single layer of IL on the surface at a relatively high temperature (~400 K). The peak height decreased with increasing temperature as can be observed in Fig. 13, which indicates that the ordering of the molecules near the surface decreased with increasing temperature. In addition to the strong layering at the proximity of the surface, we observed weak layering at the vaporliquid interface, which was pronounced at a relatively low temperature due to the orientational ordering of cations at the vaporliquid interface.
3.7. Interfacial tensions
Fig. 12. Equilibrium contact angle of IL drop of 250 ion pairs at various temperatures.
The wetting of a surface strongly depends on the interfacial tensions between different phases along the three-phase contact line. In this direction, several computational and experimental studies have been performed to correlate the wetting behavior to the interfacial tensions.50,54,55 For example, Restolho et al. [56] correlated the surface tension with the contact angle of the IL droplet and found that the higher the surface tension, the higher the contact angle of the IL droplets. Moreover, Gao et al. [57] reported a similar trend. Bordes et al. [58] recently investigated the liquid surface tension values of sixteen ILs and the contact angles of IL droplets on a graphite surface and found that the affinity of the IL toward the graphite surface was better described by the solid-liquid interfacial energy than by the liquid surface tension. In view of that,
8
S. Bhattacharjee, S. Khan / Fluid Phase Equilibria 501 (2019) 112253
Fig. 13. Density profile along the Z-direction for the different temperatures (a) 250 ion pairs of [DMIM][BF4] and (b) 250 ion pairs of [PrMIM][BF4].
we calculated the surface tension as well as the solid-liquid interfacial tension [34] for both ILs at 300 K, as given in Table 1. The value of the surface tension obtained from the simulation was very close to the experimental results, as can be observed for [PrMIM][BF4]. We observed that as the chain length increased from [C1eC3MIM] [BF4], the contact angle decreased with the value of gLV , which is in line with the work of Bordes et al. [58] Interestingly, the solid-liquid interfacial tension as well as the work of adhesion decreased as the chain length increased. Therefore, the decrease in the contact angle as the chain length increased was also owing to the graphite-IL interfacial tension. 3.8. Cylindrical drop Generally, there is a significant system size effect on spherical droplet due to the finite length of the three-phase contact line at a nano-scale; thus, the line tension plays an important role in the contact angle of the droplet. Therefore, modified Young's equation is used to determine the macroscopic contact angle of the droplet by extrapolating the microscopic contact angle of droplets of various sizes [59-64]. However, this method is very expensive for viscous liquid, such as ILs. Moreover, the base layer of the IL droplet was found to be non-spherical for a spherical drop, as shown in the snapshot of the 1st layer near the surface (Fig. 9), due to the heterogeneity in the structure. Thus, it is very difficult to obtain an accurate contact angle for the spherical droplet. Moreover, the irregular shape of the three-phase contact line may affect the ordering of the molecules near the surface, which may lead to improper wetting behavior of the ILs if the droplet size is relatively small. Conversely, the system size effect in the 2D cylindrical drop was lesser than that in the spherical droplet, as can be observed in several works [61,63]. For example, Scocchi et al. [61] compared the contact angle of the water droplet on the graphite surface for the spherical and cylindrical filament, in detail. They found that the contact angle obtained from the cylindrical filament matched well with the macroscopic contact angle of the spherical droplet obtained from the modified Young's equation. Additionally, they observed that the modified Young's equation is not applicable for small-sized droplets. Therefore, expectedly, obtaining the macroscopic contact angle through the modified Young's equation would
be very difficult for ILs. Further, Weijs et al. [63] studied the effect of line tension for LennardeJones (LJ) nano-droplet by comparing 3D spherical shaped and 2D cylindrical shaped droplets and found that the line tension does not affect the cylindrical droplet. A similar observation has also been found in several works [60-62]. However, these works are limited to only simple fluids. In this study, a cylindrical filament of IL was placed on the surface to create a onedimensional droplet in the xey plane as infinite. The contact angle of the filament was determined by locating the vapor-liquid interface in the other direction and comparing with the contact angle obtained from the spherical droplet. Fig. 14a represents the typical snapshot of the equilibrium cylindrical drop of [PrMIM][BF4] on the surface, where the y-dimension is infinite, and Fig. 14b shows the snapshot of the 1st layer of [PrMIM][BF4] near the
Fig. 14. (a) Cylindrical filament of IL [PrMIM][BF4]. (b) Snapshot of first layer [PrMIM] [BF4].
Table 1 The results of surface tension and work of adhesion of ILs at temperature 300 K. Ionic Liquid
gLV (mN/m)
Contact angle (q )
gSL (mN/m)
Literature, gLV (mN/m)
Work of adhesion Wad ¼ gLV ð1 þcos qÞ (mN/m)
Experimental [DMIM][BF4] [PrMIM][BF4]
57:7±3:1 49:5±4:2
71±1 62±1
35:9±1:0 31:5±1:9
e 47:0 [53]
77:2±5:5 72:5±5:8
S. Bhattacharjee, S. Khan / Fluid Phase Equilibria 501 (2019) 112253
surface. The same system can also be studied in a different way, in which the x-dimension would be infinite and the contact angle would be evaluated from the other dimension. These two systems can be used to determine the contact angle of the heterogeneous system, in which the contact angle of the droplet in the x- and ydirections are significantly different, as can be observed from the previous studies [19,21]. However, no significant difference was observed for the contact angle in the x and y dimensions for the studied ILs. However, quite expectedly, the contact angle for a more heterogeneous system such as [BMIM]/[HMIM] on polar surfaces (TiO2/SiO2) would be different, which is a topic for future work. Fig. 15a shows that the contact angle obtained from the spherical drop as well as from the cylindrical drop for droplets of different sizes were found to be almost identical when the drop consisted of more than 200 ion-pairs of [DMIM][BF4]. For small spherical droplets (ion pair < 200), the contact angle continuously increased with the droplet size. Conversely, the contact angle of the cylindrical droplet remained almost unchanged even for relatively small droplets (Fig. 15). Therefore, the cylindrical droplet could be useful in determining the wetting behavior for bulky and asymmetric molecules without much difference in contact angles compared to the case with the spherical droplet. However, the contact angle for [PrMIM][BF4] in the cylindrical drop was found to be higher (5e6 ) than that of the spherical drop, as shown in Fig. 15a. This difference in contact angle may be due to the system size effect (250 ion pair may not be sufficient for the spherical droplet of [PrMIM][BF4]), the irregular shape of the three-phase contact line, etc. Overall, the contact angle evaluated from the cylindrical drop was very close to that of a spherical droplet. Both the methods suggest that the contact angle of the [PrMIM][BF4] droplet was lesser than that of [DMIM][BF4]. Therefore, the length of the alkyl chain plays an important role in the wetting behavior of the ILs. The density profile for the cylindrical droplet was also evaluated and compared with that of the spherical droplet, as shown in Fig. 15b. It can be observed that the density profiles of [DMIM][BF4] obtained from both methods are very similar. Similarly, the orientation profiles obtained from both methods are very similar (not shown). However, it can be observed from the snapshot of the 1st layer of the IL molecules (Figs. 9 and 14) on the surface that the direction of the alkyl chain of cations are normal to the interface at the three-phase contact line, which consequently affects the arrangement of the IL molecules in subsequent layers inside the three-phase contact line. The curvature of the three-phase contact line was quite different for the spherical droplet compared to that of the cylindrical droplet. Therefore, quite expectedly, the packing of molecules would be different inside the three-phase contact line for both methods; therefore, the equilibrium contact angle of the
9
cylindrical droplet would be affected. However, the system size effect for the cylindrical droplet was very negligible compared to that for the spherical droplet, in accordance with the literature [6063]. 4. Conclusions The wetting behaviors of two different ILs with imidazoliumbased cation were investigated using MD simulation to understand the role of the alkyl chain attached to the imidazolium ring of the cations. The density profile and orientation order parameter profile for a spherical droplet of both the ILs showed the strong layering of cations and anions near the surface, which is also in line with other ILs [19e21]. The alkyl groups of [PrMIM] in the vaporliquid interface as well as in the three-phase contact line were found to be oriented towards the vapor phase. The contact angle of the droplet was found to decrease with increasing temperature for both ILs, resulting in a thin film on the surface at a relatively high temperature. The contact angle of the [PrMIM[BF4] droplet was found to be lesser than that of [DMIM][BF4] for the employed temperature range. The behavior of the contact angle of the IL droplet in response to increases in temperature and alkyl chain length of cations strongly depends on the surface tension and the graphite-IL interfacial tension. The surface tension and graphiteeIL interfacial tension decreases as the length of the alkyl chain increases, thereby enhancing the wetting properties of the ILs with relatively long alkyl chains. However, our studies considered only two ILs with short alkyl chains (C1eC3MIM) of cation on the graphite (non-polar) surface. As both surface tension and graphiteeIL interfacial tension decreases with increasing alkyl chain length, it is quite expected that longer alkyl chain (>C3MIM) of the cation would behave similarly, which is also the observed case for imidazolium-based cation with a different anion (NTF2) on the graphite surface [58], PTFE] [1], etc. The solideIL interfacial tension also plays an important role in the wetting behavior of IL; consequently, the effect of the alkyl chain can be different for other surfaces [1] (particularly on polar surfaces), and further investigation is required, in detail. The wetting behaviors of both ILs were also investigated for cylindrical droplet and compared with the contact angle obtained for the spherical droplet. The contact angles of [DMIM][BF4] were found to be almost identical for the spherical and cylindrical droplets for large-sized droplets (>200 ion pair for the studied ILs). However, for relatively small droplets, the contact angle of the spherical droplet continuously decreases with a decrease in the systems size. Conversely, the contact angle obtained from the cylindrical droplet was approximately the same for all the system sizes. Therefore, the cylindrical droplet was found to be efficient for
Fig. 15. (a) Contact angle of [DMIM][BF4] and [PrMIM][BF4] for spherical vs cylindrical droplet. (b) Compared Density profile along the z-axis for cylindrical and spherical droplet for 250 ion pairs, here dotted lines representing cylindrical filament and solid lines for the spherical filament.
10
S. Bhattacharjee, S. Khan / Fluid Phase Equilibria 501 (2019) 112253
characterizing the wetting behaviors of asymmetric and bulky molecules at the nanoscale considering the anisotropic behaviors of ILs in the x and y dimensions as well as the system size effect. Acknowledgments We would like to acknowledge the Department of Chemical and Biochemical at IIT Patna for providing computational resources for this research. This work is supported by grants from Science and Engineering Research Board (SERB), Government of India (Project File no. ECR/20l8/002600). We also would like to acknowledge the Center for Development of Advanced Computing (C-DAC), India for generous allocation of computing resources. We thank Dr. Snehasis Daschakraborty and Dr. Atanu Metya for a very fruitful discussion and suggestions.
[23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39]
References [1] M.M. Pereira, K.A. Kurnia, F.L. Sousa, N.J.O. Silva, J.A. Lopes-da-Silva, J.A.P. Coutinho, M.G. Freire, Phys. Chem. Chem. Phys. 17 (2015) 31653. [2] M. Poleski, J. Łuczak, R. Aranowski, C. Jungnickel, Physicochem. Probl. Miner. Process 49 (2013) 277. [3] E.J. Maginn, Acc. Chem. Res. 40 (2007) 1200. [4] A. Noda, K. Hayamizu, M. Watanabe, J. Phys. Chem. B 105 (2001) 4603. [5] H. Tokuda, S. Tsuzuki, M.A. Susan, K. Hayamizu, M. Watanabe, J. Phys. Chem. B 110 (2006) 19593. [6] T. Welton, Chem. Rev. 99 (1999) 2071. [7] J.L. Anderson, J.K. Dixon, J.F. Brennecke, Acc. Chem. Res. 40 (2007) 1208. [8] M.H. Kowsari, S. Alavi, M. Ashrafizaadeh, B. Najafi, J. Chem. Phys. 129 (2008) 224508. [9] S. Nagelberg, L.D. Zarzar, N. Nicolas, K. Subramanian, J.A. Kalow, V. Sresht, D. Blankschtein, G. Barbastathis, M. Kreysing, T.M. Swager, M. Kolle, Nat. Commun. 8 (2017) 14673. [10] N.V. Plechkova, K.R. Seddon, Chem. Soc. Rev. 37 (2008) 123. [11] A. Somers, P. Howlett, D. MacFarlane, M. Forsyth, Lubricants 1 (2013) 3. [12] T. Batchelor, J. Cunder, A.Y. Fadeev, J. Colloid Interface Sci. 330 (2009) 415. [13] R. Burt, G. Birkett, M. Salanne, X.S. Zhao, J. Phys. Chem. C 120 (2016) 15244. [14] G.V.S.M. Carrera, C.A.M. Afonso, L.C. Branco, J. Chem. Eng. Data 55 (2010) 609. n, T.J. Wynn, Z.M. Marcin, J. Phys. Chem. B 118 (2014) 3661. [15] H.J. Castejo [16] I. Delcheva, D.A. Beattie, J. Ralston, M. Krasowska, Phys. Chem. Chem. Phys. 20 (2018) 2084. [17] I. Delcheva, J. Ralston, D.A. Beattie, M. Krasowska, Adv. Colloid Interface Sci. 222 (2015) 162. [18] X. Gong, B. Wang, A. Kozbial, L. Li, Langmuir 34 (2018) 12167. [19] Y. Guan, Q. Shao, W. Chen, S. Liu, X. Zhang, Y. Deng, J. Phys. Chem. C 121 (2017) 23716. [20] C. Herrera, G. García, M. Atilhan, S. Aparicio, J. Phys. Chem. C 119 (2015) 24529. [21] S. Malali, M. Foroutan, J. Phys. Chem. C 121 (2017) 11226. [22] S. Millefiorini, A.H. Tkaczyk, R. Sedev, J. Efthimiadis, J. Ralston, J. Am. Chem. Soc. 128 (2006) 3098.
[40] [41] [42] [43] [44] [45] [46] [47] [48] [49] [50] [51] [52] [53] [54] [55] [56] [57] [58] [59] [60] [61] [62] [63] [64]
M. Paneru, C. Priest, R. Sedev, J. Ralston, J. Phys. Chem. C 114 (2010) 8383. A. Quinn, R. Sedev, J. Ralston, J. Phys. Chem. B 107 (2003) 1163. K.S. Rane, J.R. Errington, J. Chem. Phys. 141 (2014) 174706. J. Restolho, J.L. Mata, B. Saramago, J. Phys. Chem. C 113 (2009) 9321. F. Taherian, F. Leroy, L.-O. Heim, E. Bonaccurso, N.F.A. van der Vegt, Langmuir 32 (2016) 140. W.-K. Kuo, J.-J. Hsu, C.-K. Nien, H.H. Yu, ACS Appl. Mater. Interfaces 8 (2016) 32021. F. Ribet, L. De Pietro, N. Roxhed, G. Stemme, Sens. Actuators B Chem. 267 (2018) 647. S. Chen, J. Wang, T. Ma, D. Chen, J. Chem. Phys. 140 (2014) 114704. S. Khan, J.K. Singh, Mol. Simul. 40 (2013) 458. A. Kozbial, C. Trouba, H. Liu, L. Li, Langmuir 33 (2017) 959. M. Lundgren, N.L. Allan, T. Cosgrove, N. George, Langmuir 18 (2002) 10462. A.K. Metya, S. Khan, J.K. Singh, J. Phys. Chem. C 118 (2014) 4113. M. Paneru, C. Priest, J. Ralston, R. Sedev, J. Adhes. Sci. Technol. 26 (2012) 2047. dua, J. Phys. Chem. B 108 (2004) 16893. J.N. Canongia Lopes, A.A.H. Pa dua, J. Phys. Chem. B 110 (2006) 19586. J.N. Canongia Lopes, A.A.H. Pa W.L. Jorgensen, D.S. Maxwell, J. Tirado-Rives, J. Am. Chem. Soc. 118 (1996) 11225. S. Maolin, Z. Fuchun, W. Guozhong, F. Haiping, W. Chunlei, C. Shimou, Z. Yi, H. Jun, J. Chem. Phys. 128 (2008) 134504. Q. Dou, M. Sha, H. Fu, G. Wu, ChemPhysChem 11 (2010) 2438. M. Sha, G. Wu, H. Fang, G. Zhu, Y. Liu, J. Phys. Chem. C 112 (2008) 18584. M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Clarendon Press, 1987. L. Martínez, R. Andrade, E.G. Birgin, J.M. Martínez, J. Comput. Chem. 30 (2009) 2157. S. Plimpton, J. Comput. Phys. 117 (1995) 1. J.H. Walther, J. Comput. Phys. 184 (2003) 670. D.J. Evans, B.L. Holian, J. Chem. Phys. 83 (1985) 4069. W. Humphrey, A. Dalke, K. Schulten, J. Mol. Graph. 14 (1996) 33. M.J. de Ruijter, T.D. Blake, J. De Coninck, Langmuir 15 (1999) 7836. J.Z. Yang, X. Wu, X. Li, J. Chem. Phys. 137 (2012) 134104. A. Ghoufi, P. Malfreyt, D.J. Tildesley, Chem. Soc. Rev. 45 (2016) 1387. S. Wang, Y. Zhang, N. Abidi, L. Cabrales, Langmuir 25 (2009) 11078. B.L. Bhargava, S. Balasubramanian, J. Am. Chem. Soc. 128 (2006) 10073. W.-G. Xu, L. Li, X.-X. Ma, J. Wei, W.-B. Duan, W. Guan, J.-Z. Yang, J. Chem. Eng. Data 57 (2012) 2177. G. Vakili-Nezhaad, M. Vatani, M. Asghari, I. Ashour, J. Chem. Thermodyn. 54 (2012) 148. P.J. Carvalho, M.G. Freire, I.M. Marrucho, A.J. Queimada, J.A.P. Coutinho, J. Chem. Eng. Data 53 (2008) 1346. J. Restolho, J.L. Mata, B. Saramago, J. Colloid Interface Sci. 340 (2009) 82. L. Gao, T.J. McCarthy, J. Am. Chem. Soc. 129 (2007) 3804. E. Bordes, L. Douce, E.L. Quitevis, A.A.H. Padua, M. Costa Gomes, J. Chem. Phys. 148 (2018) 193840. R.C. Dutta, S. Khan, J.K. Singh, Fluid Phase Equilib. 302 (2011) 310. H. Peng, G.R. Birkett, A.V. Nguyen, Mol. Simul. 40 (2013) 934. G. Scocchi, D. Sergi, C. D'Angelo, A. Ortona, Phys. Rev. E. 84 (2011), 061602. D. Vanzo, D. Bratko, A. Luzar, J. Chem. Phys. 137 (2012), 034707. J.H. Weijs, A. Marchand, B. Andreotti, D. Lohse, J.H. Snoeijer, Phys. Fluids 23 (2011) 022001. T. Werder, J.H. Walther, R.L. Jaffe, T. Halicioglu, P. Koumoutsakos, J. Phys. Chem. B 107 (2003) 1345.