Wetting behaviour of a translating sessile nanodrop under electrostatic actuation

Abstract

Molecular dynamics (MD) simulation is performed to investigate the wetting characteristics of a nanosized pure water droplet subjected to a differential electric field. The drop considered is placed on strips of silicon electrode which can be charged and switched progressively in a direction. The results of the switching electrodes show that the droplet translates over the substrate and its dynamics can be controlled by tuning the electrode actuation. The wetting phenomenon during translation shows two distinct stages in which a precursor film forms first which subsequently drags the bulk liquid with the help of progressive switching. The frequency of the shifting charged region below the drop and magnitude of the assigned charge both have a significant impact on the translation and can be optimized for desired translation parameters. Controlled mobility of the nanodrops of a polar liquid using an electric field is perfectly aligned with the rapid technological development in nano-mission and may open up applications in the areas of biomedical and applied chemistry research.

---

1. Introduction

Dynamic wetting under the influence of an electric field or electrowetting has attracted considerable research interest for its promising aspects in small-scale liquid handling. After its discovery by Lippman¹ and substantiation by Berge,² electrowetting has been used in a vast area of engineering and biomedical applications like hot spot cooling of electronic chips,³ liquid lenses,⁴ electro wetting displays,⁵ biomedical micro-mixtures⁶ etc. The macroscopic nature of wetting under the electric field is modelled by well-established formulation proven by a vast range of experimental^7–15 and theoretical studies^16–22 and practiced extensively in various micro scale applications.^1–6 However, the growing demand of miniaturization and the rapid advancement of nanotechnology raise the research attention towards nano-scale wetting behaviour under the influence of external electric field. At the nano-scale the presence of electric field alters the circumstances of wetting from that in the macroscopic scale. As the size of the droplet is comparable to the Debye screening length, the electric field permeates more into the bulk of the droplet. As a result, the electric energy E·D_i (E is the electric field and D_i is dipole moment) has also to be accounted with the thermal energy k_BT (T is absolute temperature and k_B is Boltzmann constant) inside the droplet.^23,24 Moreover, the polarization of the surface molecules is also having an effect on surface tension and contact angle.^23,25,26

MD studies have been carried out in order to explore the microscopic wetting phenomenon under electric field. Daub et al.²³ studied the wetting of a pure water nano-droplet on a graphite surface under the influence of electric field by MD simulation. Their study reveals that the contact angle and the shape of the droplet depends on both polarity and the direction of the electric field as the alignment of the dipolar water molecules changes with it, causing differences in the hydrogen bond distribution at the interfacial regions. The MD study of Yuan and Zhao²⁴ illustrates the transport properties of the precursor film of a pure water droplet on a gold surface under the influence of electric field. They reported the precursor film to have a solid like behaviour and its propagation obeys power law with time. They established the precursor film to be a new length scale to address the stress singularity at the three phase contact line and developed molecular kinetic theory to explore the dynamics of moving contact line under external electric potential. Deformation of a nano-sized pure water droplet on silicon surface due to electric field parallel to the solid–liquid interface is also investigated by Song et al.²⁷ They observed that the asymmetry produced in the droplet due to electric field first increases and then decreases with the increasing field strength. The influence of surface roughness^25,28 and ion-concentration^29,30 on electro-wetting behaviour of a nano-droplet is also investigated by MD studies. A comprehensive overview of the MD studies on wetting under the electric field is discussed in the literature by Daub et al.³¹ All these previous literatures successfully captured and illustrated the microscopic mechanism of the wetting phenomenon of a nano-droplet under the influence of electric field.

However, the existing analyses were focused on the uniform electric field throughout the domain; the effect of differential electric field on the wetting behaviour of a nano-droplet is still unexplored. Moreover, open question remains whether the translation behaviour of droplets (due to differential electric field) observed in the macroscopic scale still persists in the nanoscale and how its physics differs from the larger scale. In the current research, MD simulations are performed to investigate the wetting behaviour of a sessile water nano-droplet on silicon surface under electric field gradient. The methodology of the molecular dynamic simulation is explained in the next section. Based on the simulation outcome the wetting phenomenon is analysed in the third section, before mentioning the salient conclusion.

2. Methodology

In the present study open source MD code LAMMPS³² is used to carry out the simulation of the dynamics of a nano-sized water droplet under differential electric field. An extended simple point charge (SPC/E)³³ water model is used to simulate the droplet. The charge at the oxygen and hydrogen sites are q_O = −0.8476e and q_H = +0.4238e. The O–H bond length and H–O–H angle is considered as 1 Å and 109.47° and constrained by SHAKE algorithm.³⁴ The solid molecules are modelled as diamond cubic silicon structure with a lattice constant of 5.43 Å. The interactions between the molecules are modelled by Lennard-Jones³³ and columbic forces. The potential energy due to the interaction forces between any two atoms can be expressed as

The first term in eqn (1) represents the pair wise dispersion and repulsion forces between molecules i and j at r_ij distance apart. Lennard-Jones energy (ε) and distance (σ) parameters for oxygen and hydrogen atoms are considered to be ε_O–O = 0.1556 kcal mol⁻¹, σ_O–O = 3.166 Å and ε_H–H = 0.0 kcal mol⁻¹, σ_H–H = 3.166 Å. Lorentz–Berthelot mixing rule ³⁵ is applied to determine the LJ interaction parameters between silicon and oxygen atoms. The distance and energy parameters of the solid atoms are σ_s–s = 3.408 Å and ε_s–s = 0.584 kcal mol⁻¹. Initially the charges at the solid atoms are considered to be zero. The electrostatic interactions between atoms are expressed by last term in eqn (1). q_i and q_j are the charges at the atom sites, ε₀ is the permittivity of vacuum.

The molecular dynamic simulation is carried out at a constant temperature of 300 K. Nose–Hoover thermostat³⁶ is utilized to relax the system temperature during 100 time steps. Time integration is carried out by Verlet algorithm³⁷ with a time step of 1 fs. To achieve accurate result the truncations of LJ and columbic interaction is considered more than 2.5σ.^38,39 Truncations of the LJ interactions between molecules are done smoothly at 14 Å.³⁰ We incorporate a particle–particle particle mesh solver (pppm)⁴⁰ to imitate long range electrostatic interactions. The real space cutoff for columbic interactions is kept 12.0 Å. Repeated check shows that the increment of the cut-off radius does not have much effect beyond the above mentioned values. Periodic boundary conditions are implemented in all the boundaries. The size of the simulation domain in the z axis is kept high (1000 Å) in order to minimize the interactions between the periodic images along it.

The simulation domain contains 3000 water molecules and 10 000 solid molecules. In order to attain a shape of a droplet from the water molecules a separate simulation is carried out, where water molecules are stacked in a cubic structure on silicon like surface and equilibrated in NVT ensemble for 200 000 time steps. Then it is placed in the main simulation domain over a silicon surface and equilibrated for another 100 000 time steps before the application of charge in the substrate. To reduce computational cost the solid wall is considered to be frozen.

The method of Liu and Xu⁴¹ is adopted to determine the contact angle of the moving droplet. However, an elliptic fit (instead of circular fit) to the outermost molecule positions (at the vertical cross section) is used to obtain the dynamic contact angles at front and rear end of the moving droplet owing to the asymmetry produced in the drop shape by non-uniform electric field. Due to the presence of charge the shape of the droplet is departed from spheroidal nature at the vicinity of the solid surface. Therefore, the contact angle is determined by elliptic fit to the spheroidal portion only (neglecting the incompliant region).²³

3. Results and discussion

In order to mimic the macroscopic electrowetting setup for droplet translation the solid wall is divided in to several electrically disconnected strips like regions with a span equal to 10.86 Å as shown in Fig. 1. The charge is assigned to two consecutive strips at a time at one end of the droplet to generate a deferential electric field. After a predefined time interval the electrification of the strip is shifted in a manner shown in the Fig. 1 as cycle 1 and 2 respectively. In the first cycle strip 1 and 2 is charged; after a specified time period (in cycle 2) the electrification of strip 1 is switched off and strip 3 is given charge. The process is repeated till the last strip. In the present study the charge switching has been done in one time step (∼fs). This could be a deviation from the potential realistic execution where a time lag would be associated with the switching mechanism. However, the study provides an insight to the behaviour of a translating nano-droplet actuated by external electric field in ideal situation. The macroscopic experiments^7–9 and the simulations^16–18 based on continuum hypothesis show a displacement of the mass centre of the droplet when subjected to a differential electric field. A similar behaviour of translation of a nano-sized droplet is observed in the present molecular dynamic study.

Fig. 1 Snap shot of the simulation domain. The green and blue strips represent the individual electrically disconnected divisions of the solid surface where charge is assigned during a specified time period (cycle). Two consecutive strips are electrified at a particular cycle and shifted with time by keeping an overlapped strip with the previous cycle.

Fig. 2 represents the droplet locations at different time step during its motion as a result of five consecutive electrode switching. From the Fig. 2 it is clearly evident that the contact angle changes at both front and rear end of the droplet during its movement due to switching of charged region. The variation of contact angles is plotted as a function of time in the Fig. 3 for the first two cycles (0 ns to 0.8 ns). The time period (τ) of the electrification in each cycle is 400 ps and the charge concentration on the solid surface is +0.008e. Initially, in the absence charge both the leading and trailing end of the droplet show approximately same contact angles (leading edge contact angle = 86.99° and trailing edge contact angle = 87.1°) due to the statistically symmetric random motion of the molecules. The contact angle of pure water on silicon surface in the present simulation is in a close agreement with the earlier reported values of 86.47°.²⁷ With the application of electric field the contact angles at both the sides of the droplet reduces in accordance with the previous microscopic studies.^23,28 However, at the rear end of the droplet the influence of electric field is lesser compared to the front end as the charge is assigned locally at the front end of the droplet in order to pertain a differential electric field. In each cycle initially the contact angle reduces very fast and then it increases with time before the switching of electrode. As the droplet is translating during each cycle, the variation of the contact angle is plotted as a function of the displacement (of the mass center of the droplet) in Fig. 4 for two consecutive actuation cycles (0 ns to 0.8 ns). It can be observed from Fig. 4 that, during the decrement of the contact angle the droplet travelled a finite distance. The displacement is higher during the growth of the contact angle.

Fig. 2 Side view of the droplet at different time instants (t = 0.05 ns, 0.5 ns, 0.95 ns, 1.4 ns and 2.0 ns) during translation. The charged region at each time instant is shown by a brace at the bottom of the solid surface. Charge concentration and time period are kept representatively as +0.008e and 400 ps.

Fig. 3 Temporal evaluation of contact angle during the first two cycles (t = 0.0 ns to 0.8 ns). The contact angles are measured at a time interval of 0.02 ns. The points are joined by dotted line for the guidance of eye. The switching of the charged region is shown by black dashed lines. The insets represent the snapshots of the droplet at the instants pointed by the arrows.

Fig. 4 Variation of contact angle as a function of displacement during the first two cycles (t = 0.0 ns to 0.8 ns). The contact angles are measured at a time interval of 0.02 ns. The clustering of data points are marked by blue dashed ellipses and distinguished as region A, B and C. The direction of the electric field is schematically shown at the inset.

When charge is assigned at the solid surface the dipolar water molecules try to reorient along the electric field by keeping the H atom in the direction of the electric field and negatively charged O atom against it. As positive charge is applied below the front end of the droplet the electric field constrained the molecular orientation with the H atom pointing outwards from the solid surface facing the bulk fluid. This causes an increment in the hydrogen bonds near the interfacial region at the front end of the droplet. But this hydrogen bond once again neutralizes as the bulk motion of the drop happens. As a result one can see the oscillatory nature of hydrogen bond counting.

Fig. 5 represents the variation of the number of hydrogen bond with time during the charge assignment at the solid surface. The number of hydrogen bond has a strong influence on interfacial tension; an increment in the hydrogen bond lowers the interfacial tension.²³ This causes a reduction of the contact angle at front end of the droplet.

Fig. 5 Temporal evaluation of the number of hydrogen bond before and after the assignment of charge at the solid surface. The insets represent the close view of the structure of the droplet front end at A and B. The blue dotted lines in the inset represent the hydrogen bond formed between water molecules.

To dig down further we analyzed the motion of the water molecules (in response to the electric field) at different region of the droplet by tracking their trajectories. Fig. 6 shows the representative of the path lines of the molecules at different regions in the droplet during time interval of 0.4 ns 0.8 ns (cycle 2). The initial positions of the molecules are shown by a spherical marker. It is observed that the path lines of molecules near the solid surface congregated at small region. This signifies low diffusion of molecules due the constraint imposed by the solid–liquid interaction. However, the unimpeded water molecules at the liquid–gas interface dispersed quickly towards the actuated region due to electric field gradient producing longer path lines before reaching to the solid surface. This observation is in agreement with the result reported by Yuan and Zhao²⁴ regarding propagation of precursor film under electric field. The rapid dispersion of the water molecules lowers the contact angle by a significant amount at the front end of the droplet. The difference in the contact angles between the front and rear end of the droplet perturbs its equilibrium and a translation of the centre of mass of the droplet occurs to regain its minimum energy condition. The difference between the contact angles at the front and rear end reduces with translation. In summary, the translation of the droplet occurs in two stages: there is an initial displacement due to rapid dispersion of the molecules from the liquid–gas interface at the immediate vicinity of the charged region; followed by the displacement of mass owing to the difference of contact angle between front and rear end. The above discussion indicates that, the translation behaviour of the droplet due to external electric field does not obey a single law. As it shows two distinct behaviours during the displacement, mathematical modelling of the phenomenon can confirm the observed fact and provides further insight to physics involved with it. First stage can be explained mathematically through molecular kinetic theory proposed by Erying et al.⁴² Here the spreading velocity can be expressed by

Fig. 6 Path lines of the representative molecules from different regions of the droplet during its translation in a complete cycle (t = 0.4 ns to 0.8 ns). The coordinate along vertical direction is given in bigger scale for clarity. The spherical marker in the graph represents the coordinates at the beginning of the cycle. The initial positions of the molecules with respect to the droplet are shown in the inset with corresponding colors.

In this expression, κ₀ is the frequency of the molecular jump between surface sites at a distance λ. k_B is the Boltzmann constant and T is absolute temperature. n is the density of sites on the silicon surface. The driving work per unit area is expressed as w. The expression in eqn (2) can further be modified following the work of Derjaguin and Churaev⁴³ and with the inclusion of average electrical energy due to the interaction of the dipole and electric field^24,44

where w_V work done per unit area by the Van der Waals interactions, w_S is the work done by structural force, w_P is the work per unit area due to polar interaction between water molecules and w_E is average electric energy due to the presence of electric field, which can be expressed as^24,44where L(x)(= coth(3x) − 1/(3x)) is the Langevin function. The quasi steady velocity of the second stage of the translation can be obtained by balancing the driving and opposing forces during the motion of the droplet. The driving force for droplet motion due to the presence of electric field can be expressed ashere, ds is the elemental length along the circumference (C) of the droplet foot print. n⃑ and i⃑ are the unit vector normal to ds and along the direction of motion respectively. γ and θ are the liquid–gas interfacial tension and dynamic contact angle. q_j is the charge at the atom sites. The opposing forces are mainly composed of the hydrodynamic drag by the solid surface and the contact line friction force. Subramanian et al.⁴⁵ determined the hydrodynamic resistance to a spherical cap shaped droplet (provided by solid surface) based on lubrication theory as

F_R = −6πμv_xR[g(θ, 1 − ε) − g(θ, 0)] (6)

where, . ε is the ratio of slip-length and droplet radius (R). v_x and μ are the velocity at the direction of motion and viscosity. The energy dissipation at the contact line can be expressed as^20,46,47

F_c = ςCv_x (7)

where, ς is the coefficient of contact-line friction. From eqn (5)–(7) the governing equation of the second stage of the droplet motion can be obtained as

As the contact angle in Fig. 4 is plotted at regular time interval (0.02 ns), a cluster of data points (can be observed at region A, B and C) signifies the pinning of the droplet. Since, the charge is applied to a finite zone (at two consecutive strips) on the silicon surface; the electric field changes its direction at half the span of the charged region. The schematic representation of the electric field is shown at the inset of Fig. 4 (however, the presence of the droplet will perturb the symmetric structure of the electric field). The change in direction in the electric field restricts the molecules to move further and causes pinning of the droplet (at region A and B). A similar behaviour is also observed in macroscopic electrowetting case for sessile droplet.¹¹ The molecules of the advancing front at the immediate vicinity of the solid surface got pinned to it due to higher columbic attraction and show low dispersion. These less mobile molecules created a hindrance to the forward motion and produced a clustering of data points at the region C. A similar behavior of translation is observed in the consecutive cycles of electrode switching. There is a minimum distinction between the rear and front end contact angles persists in each cycle owing to the equilibrium of the droplet with the non-uniform electrostatic field.

From the foregoing discussion it can be argued that the frequency of the electrode switching has a strong influence on droplet displacement. Fig. 7 represents the displacement of the centre of mass of the droplet as a function of time for different switching frequency of the charged region. The charge concentration in the electrified region is fixed at +0.008e, whereas the time period (τ) of the electrification of the specified region in each cycle (before switching to the next region) is varied from 100 ps to 400 ps. Five consecutive cycle of the shifting of the charged region is performed in each case. It can be observed from the figure that, the displacement of the centre of mass is maximum for the case with 400 ps time period. The displacement declines with the decreasing τ. For 400 ps time period, during each cycle the time span suffices the complete displacement of the centre of mass (before it get pinned due to the change in direction of the electric field) by the two stages mentioned earlier. However, for a cycle with a time period of 100 ps, the front end propagates rapidly with the application of charge but the electrification shifted to the next region causing the propagation of the wetting front further, before the completion of the second stage of displacement in the previous cycle. Thus the droplet gets elongated in the direction of motion. For the time period of 100 ps the length of the droplet footprint increases in steps (shown in the inset of the Fig. 7) during the switching of the charged region. Whereas, for the case with T = 400 ps, the droplet elongates at the beginning due the fast propagation of the spreading front and then the length get reduced during the translation of the bulk (at the 2nd stage). The phenomenon is repeated at each cycle. The increase in the length of the droplet footprint causes lesser displacement of the mass centre for τ = 100 ps than in the case of τ = 400 ps. Due to the complicated nature of the eqn (3) and (8), it difficult to have an analytical solution to obtain the displacement with time. The situation become more involved when there is a co-occurrence of both the stages of translation at the same time as in the case of τ = 200 ps during latter cycles. Thus the time variation of displacement is approximately correlated by a power law fit (displacement, L = a(q_c) + b(q_c) × t^c(q_c)) as shown in Fig. 7.

Fig. 7 Displacement of the droplet mass center as a function time during translation for different time period of electrification. The simulation is performed by keeping the charge concentration +0.008e at the electrified region. The solid lines represent the power law fit to the data points. Time evolution of the droplet length during translation is shown at the inset. The snapshots of the droplet at the end of 5 consecutive cycles in each case are represented at inset (I)–(III).

At this point it is worthwhile to mention the influence of charge concentration on the translation behaviour of the droplet. Fig. 8 shows the top view of the droplets actuated by different charge concentrations assigned to the solid surface at different time instants during their motion. The charge concentration is varied from +0.006e to +0.012e. The frequency of the switching of the charged region is kept constant (τ = 200 ps). As discussed earlier that the translation of the droplet is a complex interplay between the interfacial tension, electrostatic force and the interaction of the molecules with the solid surface. For the droplet actuated by +0.006e, the change in contact angle at the front end of the droplet is smaller. The dispersion of the molecules from the liquid–vapour interface is also less. As consequence the response of the droplet to the electric field is slower. Thus the droplet could not follow the fast electrode switching after showing some initial displacement. On the other hand, for +0.012e actuation the molecules from the liquid–vapor at the immediate vicinity of the charged solid surface moves rapidly towards the charged region. This results in a fast propagation of the precursor film at the front region. Moreover, due to the higher electric field more numbers of water molecules reorient themselves keeping the hydrogen atom towards the bulk, which causes an increase in the hydrogen bond near the solid surface and consequent reduction of interfacial tension. Thus a spreading is observed at the front end which also shifts at a high velocity with the charged region; the length of the droplet however gets increased due to the fast electrode switching. For the actuation with +0.008e the situation in between the above two cases and ended up with a sharp nose like structure at front end after the translation. It has to be mentioned that, a complete translation can be observed in all the above scenarios at sufficiently large time period of electrification.

Fig. 8 Top view of the droplets at different time instant actuated by different charge concentrations; (a) q_c = +0.006e, (b) q_c = +0.008e and (c) q_c = +0.012e. The simulations are performed by the shifting the charged region at a interval of 200 ps.

Fig. 9 shows the time variation of displacement for different actuation. The time period of the electrification in each cycle is considered to be 200 ps. As expected the droplet with +0.012e actuation has higher displacement with time. The spreading of the droplet fetches more molecules toward the front end which also bring up the displacement in case of +0.012e actuation. The time–displacement curve for +0.006e actuation became flatter as the droplet seizes its motion after some initial displacement. The length of the droplet increases at the beginning and reduced to its initial value as it misses the electrode switching. The higher rate of increase in the length with time for +0.012e actuation signifies the faster propagation of wetting front.

Fig. 9 Displacement of the droplet mass center with time during 5 consecutive cycles, actuated by different charge concentration. The charge region is shifted at time interval of 200 ps. The solid lines represent the power law fit to the data points. Time evolution of the droplet length during translation is shown at the inset.

4. Conclusion

Molecular dynamic simulation is performed to investigate the atomistic details of dynamic wetting of nano sized water droplet on silicon surface under the influence of differential electric field. The gradient of electric field is produced by imposing charges in specified regions of the solid substrate.

The study reveals that, there is a translation of mass center of the droplet as observed in the macroscopic electro-wetting applications. The translation occurs in two prominent stages. At the beginning the water molecules from the liquid–vapor interface dispersed towards the electrified region causing an initial displacement. Then in the second stage the bulk fluid moves owing to the driving force produced due to the contact angle variation between front and rear end. The difference in the contact angle between the two sides of the droplet depends on the electrostatic force, orientation of the dipoles (causes change in hydrogen bond structure) and the spreading of the molecules in the previous stage. To obtain a cumulative displacement the charged region is shifted with time. The translation behavior is analyzed further by varying the switching frequency (of the charged region) and charge concentration. The result shows an increase in the length of the droplet with high frequency of electrode switching. It is also observed that the spreading of the precursor film is more with higher charge concentration. In all the cases the displacements approximately follow a power law time dependence.

Acknowledgements

Financial support of this work was received from Department of Science and Technology, India (Grant No: SB/FTP/ETA-84/2013).

Notes and references

1. G. Lippmann, Ann. Phys., 1873, 149, 546–561.
2. B. Berge, C. R. Acad. Sci., Ser. II: Mec., Phys., Chim., Sci. Terre Univers, 1993, 317, 157–163.
3. P. Paik, V. K. Pamula and K. Chakrabarty, Adaptive Hot-Spot Cooling of Integrated Circuits Using Digital Microfluidics: Proceedings of IMECE, 2005, p. 81081.
4. B. Berge, Liquid lens technology: principle of electrowetting based lenses and applications to imaging, in Proc. 18th IEEE International Conference on Micro Electro Mechanical Systems MEMS, 2005, pp. 227–230.
5. R. A. Hayes and B. J. Feenstra, Nature, 2003, 425, 383–385.
6. P. Paik, V. K. Pamula, M. G. Pollack and R. B. Fair, Lab Chip, 2003, 3, 28–33.
7. S. K. Cho, H. Moon and C. J. Kim, J. Microelectromech. Syst., 2003, 12, 70–80.
8. M. G. Pollack, A. D. Shenderov and R. B. Fair, Lab Chip, 2002, 2, 96–101.
9. M. G. Pollack and R. B. Fair, Appl. Phys. Lett., 2000, 77, 1725–1726.
10. J. S. Kuo, P. S. Mihalic, I. Rodriguez and D. T. Chiu, Langmuir, 2003, 19, 250–255.
11. C. G. Cooney, C. Y. Chen, M. R. Emerling, A. Nadim and J. D. Sterling, Microfluid. Nanofluid., 2006, 2, 435–446.
12. U. C. Yi and C. J. Kim, EWOD Actuation with electrode-free cover Plate. In 13th International Conference on Solid-State Sensors, Actuators and Microsystems, Seoul, South Korea, 2005.
13. J. Lee, H. Moon, J. Fowler, T. Schoellhamer and C. J. Kim, Sens. Actuators, A, 2002, 95, 259–268.
14. H. Moon and J. Kim, Sens. Actuators, A, 2006, 130, 537–544.
15. S. K. Fan, Y. W. Hsu and C. H. Chen, Lab Chip, 2011, 11, 2500–2508.
16. J. Berthier, P. Clementz, O. Raccurt, D. Jary, P. Claustre, C. Peponnet and Y. Fouillet, Sens. Actuators, A, 2006, 127, 283–294.
17. A. Arzpeyma, S. Bhaseen, A. Dolatabadi and P. Wood-Adams, Colloids Surf., A, 2008, 323, 28–35.
18. H. Aminfar and M. Mohammadpourfard, Comput. Meth. Appl. Mech. Eng., 2009, 198, 3852–3868.
19. Z. Keshavarz-Motamed, L. Kadem and A. Dolatabadi, Microfluid. Nanofluid., 2010, 8, 47–56.
20. V. Bahadur and S. V. Garimella, J. Micromech. Microeng., 2006, 16, 1494–1503.
21. S. Datta, M. Sharma and A. K. Das, Ind. Eng. Chem. Res., 2014, 53, 6685–6693.
22. T. B. Jones, J. D. Fowler, Y. S. Chang and C. J. Kim, Langmuir, 2003, 19, 7646–7651.
23. C. D. Daub, D. Bratko, K. Leung and A. Luzar, J. Phys. Chem. C, 2007, 111, 505–509.
24. Q. Yuan and Y. P. Zhao, Phys. Rev. Lett., 2010, 104, 246101.
25. T. H. Yen, Mol. Simul., 2012, 38, 509–517.
26. D. Bratko, C. D. Daub, K. Leung and A. Luzar, J. Am. Chem. Soc., 2007, 129, 2504–2510.
27. F. H. Song, B. Q. Li and C. Liu, Langmuir, 2013, 29, 4266–4274.
28. Y. P. Zhao and Q. Yuan, Nanoscale, 2015, 7, 2561–2567.
29. F. H. Song, B. Q. Li and C. Liu, Langmuir, 2014, 30, 2394–2400.
30. C. D. Daub, D. Bratko and A. Luzar, J. Phys. Chem. C, 2011, 115, 22393–22399.
31. C. D. Daub, D. Bratko and A. Luzar, Top. Curr. Chem., 2012, 307, 155–179.
32. S. Plimpton, J. Comput. Phys., 1995, 117, 1–19.
33. H. J. C. Berendsen, J. R. Grigera and T. P. Straatsma, J. Phys. Chem., 1987, 91, 6269–6271.
34. G. Ciccotti, M. Ferrario and J. P. Ryckaert, Mol. Phys., 1982, 47, 1253–1264.
35. M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids, Oxford University Press, Oxford, U.K., 1989.
36. D. J. Evans and B. L. Holian, J. Chem. Phys., 1985, 83, 4069–4074.
37. L. Verlet, Phys. Rev., 1967, 159, 98.
38. S. D. Luca, B. D. Todd, J. S. Hansen and P. J. Daivis, J. Chem. Phys., 2013, 138, 154712.
39. J. Liu, M. Wang, S. Chen and M. O. Robbins, Phys. Rev. Lett., 2012, 108, 216101.
40. T. Darden, D. York and L. Pedersen, J. Chem. Phys., 1993, 98, 10089–10092.
41. Q. Liu and B. Xu, Langmuir, 2015, 31, 9070–9075.
42. S. Gladstone, K. Laidler and H. Eyring, The Theory of Rate Processes, McGraw-Hill Book Company, New York, 1941.
43. B. V. Derjaguin and N. V. Churaev, J. Colloid Interface Sci., 1974, 49, 249–255.
44. Y. P. Zhao and Y. Wang, Rev. Adhes. Adhes., 2013, 1, 114–174.
45. R. S. Subramanian, N. Moumen and J. B. McLaughlin, Langmuir, 2005, 21, 11844–11849.
46. T. D. Blake and J. D. Coninck, Adv. Colloid Interface Sci., 2002, 96, 21–36.
47. H. Ren, R. B. Fair, M. G. Pollack and E. J. Shaughnessy, Sens. Actuators, B, 2002, 87, 201–206.

---