Commun. Theor. Phys. 66 (2016) 467–473

Vol. 66, No. 4, October 1, 2016

Study of Friction between Liquid Crystals and Crystalline Surfaces by Molecular Dynamic Simulations∗ Yong-Wen Zhang (张永文),1 Xiao-Song Chen (陈晓松),1,2 and Wei Chen (陈卫)3,† 1

CAS Key Laboratory of Theoretical Physics, Institute of Theoretical Physics, Chinese Academy of Sciences, P.O. Box 2735, Beijing 100190, China 2 School of Physical Sciences, University of Chinese Academy of Sciences, Beijing 100049, China 3

Computer Network Information Center, Chinese Academy of Sciences, P.O. Box 349, Beijing 100190, China

(Received June 13, 2016; revised manuscript received July 8, 2016)

Abstract The lubrication characteristics of liquid crystal (LC) molecules sheared between two crystalline surfaces obtained from molecular dynamics (MD) simulations are reported in this article. We consider a coarse-grained rigid bead-necklace model of the LC molecules confined between two atomic surfaces subject to different shearing velocities. A systematic study shows that the slip length of LC lubrication changes significantly as a function of the LC-surface interaction energy, which can be well described though a theoretical curve. The slip length increases as shear rate increases at high LC-surface interaction energy. However, this trend can not be observed for low interaction energy. The orientation of the LC molecules near the surface is found to be guided by the atomics surfaces. The influence of temperature on the lubrication characteristics is also discussed in this article. PACS numbers: 81.40.Pq, 83.80.Xz

Key words: liquid crystal, friction, slip length, molecular dynamics simulation

1 Introduction Recently, using liquid crystal (LC) molecules as lubricants or lubricant additives has been shown to produce a dramatic reduction of friction to ultralow values.[1−4] We have already performed a systematic molecular dynamics (MD) study on the friction force of the confined LC monolayer under different shearing and sliding conditions by considering three different atomic structures for the confining surfaces in our previous article.[5] It has demonstrated that the frictional properties of the LC molecules depend on the competition between the effect of the pattern of surface mesh and that of the imposed sliding direction. Many studies have also focused on the influence of the lubricant-surface interaction energy on friction force. The stick-slip events and periodic breaking-reforming transitions of atomic-scale capillary water bridges can be observed for thin water films confined by hydrophilic mica sheets from MD simulations. However, only smooth sliding without stick-slip events is observed for water confined by hydrophobic graphene.[6] It was also shown in experiments[7−8] that the hydrogen bond network between hydrophilic surfaces would lead to decrease of friction with increase in shear rate. On bare mica surfaces and other hydrophobic surfaces, however, the trend is inversed. Thus, controlling the lubricant-surface interaction energy could tune the interfacial friction. In order to make fur∗ Supported

ther progress in the optimization of frictional properties of LC lubricants, in this study, we continue studying the key structural and dynamical properties of confined LC molecules under shear conditions by performing MD simulations of a simplified rigid bead-necklace model of the LC molecules.[9−11] Previous study about the boundary lubrication characteristics of this model have shown that the coarse-grained bead-necklace model capture much of the structural and dynamic properties of LC lubricant.[5] We here consider LC molecules confined by surfaces with different energy parameters of LC-surface interaction under temperature changing conditions, approximating a wide variety of materials. We will show that the attractive force between LC molecules and surfaces leads to the enhanced LC density layering away from the surfaces and partial alignment of LC molecules in the first fluid layer. And the slip length is found to change significantly as a function of LC-surface interaction energy. The article is arranged as follows. The details of the MD simulations procedure are described in Sec. 2. The simulation results are presented and discussed in details in Sec. 3. We conclude the article by a few remarks in the last section.

2 Model In this study, LC molecule is treated as a coarsegrained rigid bead-necklace model.[5,9−11] Each LC molecule consists of nine interaction sites (beads). The

by the National Natural Science Foundation of China under Grant Nos. 11504384 and 11121403 and computational resources provided by Supercomputing Center of Chinese Academy of Sciences (SCCAS) † E-mail: [email protected] c 2016 Chinese Physical Society and IOP Publishing Ltd ⃝ http://www.iopscience.iop.org/ctp http://ctp.itp.ac.cn

468

Communications in Theoretical Physics

non-bonded interactions between the beads belonging to different molecules are described by pairwise-additive Lennard–Jones (L-J) 12-6 potentials, [( σ )12 ( σ )6 ] bb bb Uij = 4ϵbb − , (1) rij rij where Uij is the interaction potential for bead i and bead j separated by a distance rij , ϵbb is the well depth of the L-J 12-6 potential, and the parameter σbb is the finite distance at which the potential is zero. The subscript b denotes bead. To be convenient, a choice for the basic units during all the simulation in this study is following: unit of length, σbb ; unit of energy, ϵbb ; unit of mass, m (the mass of one bead). From these basic units, the temperature, pressure, and time are defined, respectively, as √ 3 2 . T ∗ = T kB /ϵbb , P ∗ = P σbb /ϵbb , and τ = t∗ = t ϵbb /mσbb The fundamental units σbb , ϵbb , m and the Boltzmann constant kB are set equal to one. The center-to-center distance between the closest beads in one molecule is kept at 0.6 in terms of the reduced units, so the total length of the LC molecule is equal to 4.8.

Fig. 1 (Color online) The sketch of the simulation geometry. y axis is normal to the page. Blue LC molecules are confined by two rigid crystalline plates (black atoms). The top plate is moved with a constant driving velocity V along the x direction, and the bottom plate is fixed. The distance between the two plates is kept at 20.04.

Figure 1 shows a sketch of the simulation geometry. LC molecules are confined by two rigid crystalline surfaces. Atomic structure of the surfaces is a projection of the face-centered cubic (fcc) on the plane that is perpendicular to the (100) direction. The top and bottom surfaces have the same number of the atoms 1800, and the distance between the two surfaces is kept at 20.04. The total number of LC molecules is set to 576 in all of the MD simulations. The length of simulation box in x, y and z directions are chosen to be Lx = 30.0, Ly = 30.0, and Lz = 40.0. The periodic boundary conditions are used in all the three directions. The number density ρ of the LC molecules confined by two surfaces can be calculated by ρ = 576/(30.0 × 30.0 × 20.04) ≈ 0.0319. Each molecular bead interacts with the surface atoms via the potential as

Vol. 66

Eq. (1), the parameters ϵbs , σbs are obtained by combination rules, √ √ ϵbs = ϵbb ϵss , σbs = σbb σss , (2) where the subscript s denotes surface atom. σss is always equal to 1, and ϵss is selected in the range of 0.2 to 1.8 in this study. All the MD simulations are performed in the NVT canonical ensemble (constant number of particles, constant volume and constant temperature). We consider both equilibrium and non-equilibrium properties of the confined LC molecules. The top surface is moved with a constant driving velocity V which is varied in a range between zero to 0.95 along the x direction, and the bottom surface is always kept still. A Nos´ e-Hoover thermostat is implemented when V = 0 to control the temperature in the equilibrium state. A Langevin themostatting, which is widely used in MD simulations of sheared fluids,[5,12−13] is applied in the y direction to remove viscous heating generated in the shear flow when V > 0. As for the Langevin themostatting, the equations of motion of the i-th bead are dυi,x (t) m = Fi,x (t) , (3) dt dυi,y (t) m = Fi,y (t) − mηυi,y (t) + fi (t) , (4) dt dυi,z (t) m = Fi,z (t) , (5) dt where υi,x , υi,y and υi,z are the projection of i-th bead velocity in x, y and z direction respectively. Fi,x , Fi,y and Fi,z are the projection of the net deterministic force acting on the i-th bead in x, y and z direction. fi (t) is a δ-correlated stationary Gaussian process with zero mean, satisfying ⟨fi (t)⟩ = 0 and ⟨fi (t)fj (t′ )⟩ = 2mkB T ηδ(t − t′ )δij . η is the damping factor which is set to 0.01 in our simulation. Different types of liquid crystal phases have been observed in different temperature ranges.[14] Thus, we also consider the influence of temperature on the system. Temperatures are chosen to be T ∗ = 3.5, 4.5, and 5.5 respectively. All of the MD simulations are performed with the Largescale Atomic/Molecular Massively Parallel simulator (LAMMPS)[15] in this study. After an initial 107 time steps during which the system reaches the steady state, a typical production runs of a total of 107 time steps are carried out with the integration time step of 0.001τ .

3 Results 3.1 Fixed Surfaces The results of equilibrium MD simulations when the driving velocity is zero is reported firstly. Figure 2 shows the density profiles ρ(z)/ρ of LC molecules confined by two surfaces with different energy parameters ϵss under different temperatures. As for T ∗ = 3.5, the LC molecules maintain layered structures in z direction. The molecules

No. 4

Communications in Theoretical Physics

are divided into 20 layers corresponding to 20 peaks. The average kinetic energy of the LC molecule increases as the temperature increases. The layered structures of LC molecules are destroyed by thermodynamic fluctuation at high temperatures. As it is shown in Fig. 2(a), at T ∗ = 5.5, only the LC molecules close to the surfaces are

469

in the layered structure. Increasing the LC-surfaces interaction energy leads to more LC molecules are attracted by the surfaces. It can be observed that at a low LC-surfaces interaction, ϵss = 0.2, the height of the peaks nearest the surfaces is much smaller than those of the larger interaction energy.

Fig. 2 (Color online) Reduced density profiles of the LC molecules in z direction for temperatures T ∗ = 3.5 (red lines), 4.5 (green lines) and 5.5 (blue lines) with different energy parameters (a) ϵss = 0.2, (b) ϵss = 1.0 and (c) ϵss = 1.8. The number density ρ is equal to 0.0319. The bottom surface and the top surface are located at z = 0 and z = 20.02 respectively. The dashed line represents the position of the trough following the peak nearest the bottom surface, which is located at z = 1.52. This position is denoted by a in this article.

Fig. 3 (Color online) Variation of the adsorption capacity of bottom surface as a function of ϵss for temperatures T ∗ = 3.5 (red squares), T ∗ = 4.5 (green circles) and T ∗ = 5.5 (purple triangles).

We find that the positions of the troughs following the peaks nearest the surfaces are almost unchanged in any conditions. As for the bottom surface, those troughs are always located at z = 1.52, this position is∫ denoted by a a, and we define the adsorption capacity as 0 ρ(z)dz/ρ. It can be observed from the data reported in Fig. 3 that the adsorption capacity of bottom surface increases as the ϵss increases in different temperatures, and it is saturated when ϵss is large enough. For the LC molecules confined by two surfaces with ϵss larger than 1.0, we find that the adsorption capacity increases monotonically with decreasing temperature at a fixed value of ϵss , however, we can not observe this trend when ϵss is smaller than 1.0. Surface structures have significant effects on interface friction by affecting the distribution of angles of the adsorbed LC molecules. It has been observed in the previ-

Communications in Theoretical Physics

470

ous work that LC molecules align their long axes roughly along the surface mesh and induce a number of domains with different local orientation orders.[5] We consider here the influence of the LC-surface interaction energy on the orientation of LC molecules close to the surfaces. The distributions of | cos(θ)| and | cos(φ)| for LC monolayer

Vol. 66

closest to the bottom surface are calculated with different ϵss and temperatures. In spherical coordinate system, θ angle of the unit vector of the LC molecule measures from the fixed z direction, and the φ angle of its orthogonal projection on x-y plane measures from x direction on x-y plane.

Fig. 4 (Color online) The distributions of | cos(θ)| and | cos(φ)| for LC molecules closest to the surface with the energy parameters (a) ϵss = 0.2, (b) ϵss = 1.0, and (c) ϵss = 1.8 at different temperatures.

The distributions are plotted in Fig. 4. The values of | cos(θ)| are around 0.01, and they do not depend on ϵss and temperature. Those results indicate that θ is very close to 90◦ , the LC molecules near the surfaces almost lie in the x-y plane. However, the maximum value of P (| cos(θ)|) decreases with increasing temperature at a fixed ϵss . In contrast, this value increases as ϵss increases at a fixed temperature. In general, the effect of the temperature competes with the influence of the LC-surface interaction energy. The nearest atoms in the (100) surface are arranged along the 2D vectors in the x-y plane with the angles that are equal to ±45◦ . The fraction of LC molecules with φ = ±45◦ increases with increasing ϵss for a fixed value of temperature. The results can be understood that the larger value of ϵss makes the surfaces more attractive. However, when temperature increases for a fixed ϵss , the LC molecules get less ordered due to the thermodynamic fluctuation. 3.2 Sliding Surfaces We consider hereafter the non-equilibrium properties

of the LC molecules by the shear flow simulations. The top surface is moved with a constant nonzero driving velocity V along the x direction. The velocity profiles of the LC molecules for V = 0.2 under different temperatures are plotted in Fig. 5. For ϵss = 0.2, the slip velocity profile of LC molecules at T ∗ = 3.5 coincides with that of T ∗ = 4.5, their slopes are almost zero which indicates that adjacent layers move parallel to each other with same speeds. When the temperature is T ∗ = 5.5 for ϵss = 0.2, each layer of LC molecules moves faster than the one just below it. For a higher LC-surface interaction energy, ϵss = 1.0, the velocity profiles are nonlinear for T ∗ = 3.5 and T ∗ = 4.5, a large amount of LC molecules attach to the surfaces and move with the same velocities. However, the velocity profiles are almost linear for T ∗ = 5.5. The similar behavior can be observed for ϵss = 1.8. Lubrication layers and sticky layers close to surfaces typically occur for the shear flow simulations at high LC-surface interaction energy, and the width of sticky layer becomes smaller as temperature increases until sticky layer disappears. For T ∗ = 5.5, the layered struc-

No. 4

Communications in Theoretical Physics

ture of LC molecules is destroyed, there are only lubrication layers at the LC-surfaces interface. Our results are similar as that reported in Ref. [16], where polymer solution forms lubrication layers at weakly attract surfaces, the sticky surface layers only appear for more attractive

471

surfaces. Obviously, the temperature can also change the velocity profiles from nonlinear to linear. Figure 6 shows the results of fitting the velocity profiles for three different energy parameters ϵss and different shearing velocities V of top surface at T ∗ = 5.5.

Fig. 5 (Color online) Velocity profiles vx (z) for temperatures T ∗ = 3.5, 4.5, and 5.5 with the different energy parameters (a) ϵss = 0.2, (b) ϵss = 1.0 and ϵss = 1.8. We take 20 bins in z direction to calculate vx (z). Top surface speed V = 0.2.

Fig. 6 (Color online) Velocity profiles for the energy parameters (a) ϵss = 0.2, (b) ϵss = 1.0 and ϵss = 1.8. Top surface speeds V = 0.2 (red), V = 0.5 (green) and V = 0.8 (blue). The dashed lines are linear regression lines and their extrapolation length in z direction are the slip length b. The temperature is kept at T ∗ = 5.5.

Fig. 7 (Color online) Variation of the slip length b as a function of the parameter ϵss for different shearing velocities of top surface, V = 0.2 (orange), V = 0.5 (green) and V = 0.8 (blue); The dash lines are linear regression lines.

The extrapolation of each linear fitting in z direction corresponds to a slip length b. It turns out that, as the case in earlier studies for simple liquid system,[17−18] the slip length b = η0 /λ, where η0 and λ are the viscosity and

the friction coefficient. ∫ ∞ λ can be obtained by a Green– Kubo formula, λ = 0 ⟨Fx (t)Fx (0)⟩dt/SkB T , where Fx is the x-component of a force vector acting on the surface with area S, and the average is taken over equilibrium configurations. A simplified equation is given by λ ≈ τ0 ⟨Fx2 ⟩/SkB T , where the relaxation time τ0 ∝ σ 2 /D, σ is a characteristic length of the surface and D is the fluid diffusion coefficient. In addition, the mean squared interaction force ⟨Fx2 ⟩ ∝ (ϵ/σ 2 )2 S at the thermodynamic limit, where ϵ is the L-J simple liquid-surface energy parameter. So the dependence of the slip length b on the strength of the liquid-surface interactions for simple liquids was expressed by b = η0 /λ ∝ σ 2 η0 DkB T /ϵ2 . For LC −1 molecules, b ∝ ϵ−2 bs is predicted. We then obtain b ∝ ϵss by combining Eq. (2). Figure 7 shows the results of variation of the slip length b as a function of the parameters ϵss for three different shearing velocities V = 0.2, V = 0.5, and V = 0.8. It can be observed that the slip length increases with decreasing interaction energies. The slip length b varies according to

472

Communications in Theoretical Physics

ϵ−1 ss linearly. The value of slip length b is small with strong interaction between LC molecules and surfaces. When this interaction becomes weaker, there are larger values for the slip length b. The similar behavior can also be observed

Vol. 66

in the MD simulations of polymer solution,[16] where the slip length b increases with decreasing liquid-surface interaction energy.

Fig. 8 (Color online) Variation of the slip length b as a function of the slip velocity of top surface at different energy parameters ϵss . The red lines are linear regression fitting lines and the shadowed area is 95% confidence region.

It has been shown from the experimental results on the effect of humidity for mica surfaces,[19] that as for the water confined by hydrophilic surfaces, the friction force decreases with increasing shearing velocities. This effect is basically the result of the hydrogen bond network forming between hydrophilic surfaces. However, the trend is inversed on bare mica surfaces and other hydrophobic surfaces, because there is no hydrogen bond network formation kinetics in this process. As expected, the increase of energy parameter ϵss leads to an increase in the wettability of the surface,[20] thus complementary simulations are performed to study the dependence of the slip length on shearing velocity at different energy parameters. For simulations on simple liquid system,[21] a gradual transition in shear rate dependence of the slip length from linear to highly nonlinear is observed upon reducing the strength of wall-liquid interactions. The slip length is a linear function of shear rate at high wall-liquid interaction energies, and it is a power-law function for weak interaction energies. Figure 8 shows the results for LC molecules at T ∗ = 5.5, it can be observed that the slip length depends on the shearing rate linearly at large values of ϵss , we can also observe nonlinear behavior at ϵss = 0.2, 0.3, and 0.4, but it does not comply with the power-law, the widths of 95% confidence region for the linear regression fit are much greater than those of large values of ϵss . Friction depends on shearing velocity has been confirmed. When the surface materials are different, velocity dependence of friction will be changed. The chemical nature of the surface, which can form H-bond network,

exhibit a friction that decreases with shear velocity, but if the surface can not form such networks, this behavior is opposite.[19] So exploring influence of variation of the slip length b as a function of shearing velocity by different wallfluid interaction energies is significant. When the shearing velocity changes, we get different values of the slip length b under same energy parameter ϵss (In Fig. 7). For simulations on simple liquid system,[21] a gradual transition in shear rate dependence of the slip length from linear to highly nonlinear is observed upon reducing the strength of wall-fluid interactions. The slip length is a linear function of shear rate at high wall-fluid interaction energies, it is a power-law function for weak wall-fluid interaction energies. Our results are similar to simple liquid at high wall-fluid interaction energies, for weak wall-fluid interaction energies, we also observe nonlinear behavior but it does not comply with the power-law. Figure 8 shows this result. We add a linear regression line to each picture of ϵss . For low values of ϵss = 0.2, 0.3, and 0.4, the widths of 95% confidence region for the regression fit are large, wherefore variation of the slip length as the slip velocity increases are nonlinear and no clear trend. When the energy parameter ϵss takes larger values, we can get better linear regression lines and the slip length b increases as the shearing V increases linearly.

4 Conclusion In summary, molecular dynamics simulation was applied to study the friction between liquid crystals and crystalline surfaces for different LC-surface interaction en-

No. 4

Communications in Theoretical Physics

473

ergies, temperatures, and shearing velocities. Our results show that the LC molecules nearest the surfaces exhibit significant orientational order at high LC-surface interaction energies and low temperatures, but get less ordered when the temperature increases or interaction energies decrease. Our findings reveal that the slip length varies as a

function of the LC-surface interaction energy, which can be well described though a theoretical curve. We find that the slip length increases linearly with increase in the shearing velocity at high LC-surface interaction energies, but for weak interaction energies, no significant trend can be observed.

References

[12] P.A. Thompson and S.M. Troian, Nature (London) 389 (1997) 360.

[1] R.J. Bushby and K. Kawta, Liquid Crystals 38 (2011) 1415. [2] T. Amann and A. Kailer, Wear 271 (2011) 1701. [3] T. Amann and A. Kailer, Tribol Lett. 41 (2011) 121. [4] T. Aman and A. Kailer, Tribol Lett. 37 (2010) 343. [5] W. Chen, S. Kulju, A.S. Foster, M.J. Alawa, and L. Laurson, Phys. Rev. E 90 (2014) 012404. [6] W. Chen, A.S. Foster, M.J. Alawa, and L. Laurson, Phys. Rev. Lett. 114 (2015) 095502. [7] J. Chen, I. Ratera, J.Y. Park, and M. Salmeron, Phys. Rev. Lett. 96 (2006) 236102. [8] S. Ohmishi and A. Stewart, Langmuir 18 (2002) 6140. [9] P. Tian, D. Bedrov, G.D. Simith, and M. Glaser, J. Chem. Phys. 115 (2001) 9055. [10] P. Tian and G.D. Smith, J. Chem. Phys. 116 (2002) 9957. [11] P. Tian, D. Bedrov, G.D. Smith, M. Glaser, and J.E. Maclennan, J. Chem. Phys. 117 (2002) 9452.

[13] N.V. Priezjev, J. Chem. Phys. 136 (2012) 224702. [14] A.J. McDonald and S. Hanna, Phys. Rev. E 75 (2007) 041703. [15] S. Plimpton, J. Comp. Phys. 117 (1995) 1. [16] J. Servantie and M. Muller, Phys. Rev. Lett. 101 (2008) 026101. [17] D.M. Huang, C. Sendner, D. Horinek, R.R. Netz, and L. Bocquet, Phys. Rev. Lett. 101 (2008) 226101. [18] C. Sendner, D. Horinek, L. Becquet, and R.R. Netz, Langmuir 25 (2009) 10768. [19] J. Chen, I. Rathera, J.Y. Pack, and M. Salmeron, Phys. Rev. Lett. 96 (2006) 236102. [20] T. Werder, J.H. Walther, R.I. Jaffe, T. Halicioglu, and P. Koumoutsakos, J. Phys. Chem. 107 (2003) 1345. [21] N.V. Prizezjev, Phys. Rev. E 75 (2007) 0501605.